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ABSTRACT 

We present new observations from Z-Spec, a broadband 185—305 GHz spectrometer, of five sub- 
millimeter bright lensed sources selected from the Herschel Astrophysical Terahertz Large Area Survey 
(H-ATLAS) science demonstration phase (SDP) catalog. We construct a redshift finding algorithm 
using combinations of the signal-to-noise of all the lines falling in the Z-Spec bandpass to determine 
redshifts with high confidence, even in cases where the signal-to-noise in individual lines is low. We 
measure the dust continuum in all sources and secure CO redshifts for four out of five [z ~ 1.5 — 3). In 
one source, SDP.17, we tentatively identify two independent redshifts and a water line, confirmed at 
z=2.308. Our sources have properties characteristic of dusty starburst galaxies, with magnification- 
corrected star formation rates of 10 2-3 M yr _1 . Lower limits for the dust masses (^a few 10 s M ) 
and spatial extents (~1 kpc equivalent radius) are derived from the continuum spectral energy dis- 
tributions, corresponding to dust temperatures between 54 and 69 K. In the LTE approximation, we 
derive relatively low CO excitation temperatures (< 100 K) and optical depths (r < 1). Performing 
a non-LTE excitation analysis using RAD EX, we find that the CO lines measured by Z-Spec (from 
J= 4 — > 3 to 10 — > 9, depending on the galaxy) localize the best solutions to either a high-temperature 
/ low-density region, or a low-temperature / high-density region near the LTE solution, with the opti- 
cal depth varying accordingly. Observations of additional CO lines, CO(l— 0) in particular, are needed 
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to constrain the non-LTE models. 

Subject headings: submillimeter:galaxies — galaxies: distances and redshifts — galaxies:high-redshift — 
galaxies:ISM — line:identification 



1. INTRODUCTION 

Galaxies detected by their thermal dust emission 
at submillimctcr (submm) and millimeter (mm) wave- 
lengths (A ~ 250 — 2000 fim) comprise an important pop- 
ulation of massive systems in the early Universe that are 
thought to be undergoi ng a phase of inten se star forma- 
tion in their evolution (| Brain et~aL|[2002l) . Dust grams 
within star-forming regions in these galaxies are heated 
by incident optical and ultraviolet (UV) radiation from 
young stars and thermally re-radiate this energy at far- 
infrared (far-IR) to mm wavelengths, with the peak of 
dust emissi on occurring at ~ 60 — 200 /xm in the rest- 
frame (e.g., IDale fc Heloul2002tlHwang et al.l2Q10T) . It is 
estimated that about half of all star- formation in the Uni- 
verse is heavily obscured by dust and therefore difficult to 
identify in ev en the deepest surv eys at optical/ultraviolet 
wavelengths (|Puget et alj|1996l ). 

Observations at submm/mm wavelengths sample the 
Rayleigh- Jeans tail of the thermal d ust spectrum, whic h 
rises steeply with frequency ~ v z b (jDunne 
For observations at A > 500 /urn, the climb up this steep 
spectrum with increasing redshift roughly cancels the ef- 
fect o f cosmological dim ming with increasing distance 
(e.g., iBlain et all 12001 . meaning that galaxies with a 
fixed luminosity will have roughly the same observed flux 
density at submm/mm wavelengths for redshifts between 
1 < z < 10. This allows a distance-independent study of 
dust-obscured star-formation and galaxy evolution span- 
ning the epoch of peak star formation activity in the 
Universe (z ~ 2 - 3, e.g., lHopkinsll200l . 

Although attempts to predict the sources responsi- 
ble for the Cosmic Far-Infared Backgrou nd (CFIRB) 
have b e en made long before its detection bylPuget et all 
(fl99ll (IPartridge fc Peebles! [T967t ILow fe Tucked 11 968L 
e.g., ), the population of high-redshift and heavily 
dust-obscured galaxies (submillimete r galaxies, SMGsj 
was first revealed a decade ago (ISmail et all 119971 ; 
IBarger et all 119981: iHughes et all |1998|), and is now 
considered to prod uce most of the observed CFIRB 
(|Devlin et al.l 120091 e.g., ). Several wide-area surveys at 
850 /xm - 1.2mm have been carried out since then (e.g., 
Weifi et alll2009bHAustermann et aTll2010tlCoppin et all 



20061: IBertoldi et all 120071: IScott et al.l 120081 ). mapping 
a total of ~ 4deg 2 of sky. More recently, much 
larger area surveys have been undertaken with the 
South Pole Telescope (SPT, IVieira et all [2010h at 
A = 1.4 — 2 mm, the Balloon-bo rne Large Aperture 
Submillimeter Tel escope (BLAST, IPascale et all 120081: 
iDevlin et aLll2009D at A = 250 - 500 urn and the Her- 
schel Space Observatory (|Pilbratt et all I2010D at A = 
55 — 670 urn. Mapping a total area of ~ 200 deg 2 to 
date (IPascale et alll2008l: IDevlin et alll2009t IVieira et~all 
[2010t lEales et al.ll2010D . these surveys have uncovered a 
population of rare, and unusually bright, distant galax- 
ies. Their inferred IR luminosities and high redshifts 
are consistent with a significant fraction of these ex- 
tremely bri ght submm/mm gal axies being gravitation- 
ally lensed (|Negrello et aLll2007D . but proof requires ex- 



tensive multi- wavelength follow-up campaigns. Their ob- 
served flux densities can be magnified by factors > 10 
due to lensing by intervening foreground galaxies or 
clusters, as observed i n similarly bright systems (e.g . , 
iSwinbank etall 120101: iSolomon fe Vanden BouH 120051) . 
By targeting lensed objects, we can study the proper- 
ties of typical star forming galaxies in the early Uni- 
verse that would otherwise be inaccessible due to sen- 
sitivity limitations and source confusion. The ongoing 
Herschel- Astrophysical Tera hertz Large Area Survey (H- 
ATLAS, lEales et all I2010D in the Science Demonstra- 
tion Phase (SDP) has already covered 14.4 deg 2 out 
of the ~550 deg 2 planned, resulting in ~ 6600 sources 
(jClements et all 120101: IRigbv et all I2010D with fluxes 
measured at 250, 350, and 500 fim using the Spectral and 
Photometric Imaging R eceiver (SPIRE, i Griffin et all 
[20101 IPascale et al.ll20Tol ). and fluxes at 100 and 160 /im 
obtained with the Ph otodetector Array Camera and 
Spect rometer (PACS, iPoglitsch et~al1 l2010t Hbar et all 
120101) . Given the large areal coverage, H- ATLAS can de- 
tect the brightest (i.e. rarest) distant submm galaxies 
and is the first example where the efficient selection of 
lensed g alaxies at submm wa velengths has been demon- 
strated (jNegrello et alll2010f ). 

To understand the nature of these galaxies, in par- 
ticular whether they represent a previously undiscov- 
ered population of i ntrinsically bright sources (e.g., 
IDevriendt etaLl [2010h or are relatively normal star- 
burst galaxies lens ed by foreground structures (e.g., 
INegrello et alll200l . requires both complementary data 
at other wavelengths and measurements of their red- 
shifts. However, measuring spectroscopic redshifts for 
these sources is challenging: their positional accuracy 
from submm/mm imaging is often poor due to diffrac- 
tion limitations at these long wavelengths, and they tend 
to be highly extincted by dust, making spectroscopic 
measureme nts from optical groun d-based telescopes dif- 
ficult (e.g., IChapman et al.ir2005l ). The positional uncer- 
tainty can be overcome by finding optical/infrared coun- 
terparts, or by deep interferometric observations at ra- 
dio a nd millimeter wavelengths (e.g., iDannerbauer et all 
2002;. This not only requires large observing campaigns, 
but can also introduce selection effects in determining 
the properties of the SMG population. In particular, 
the combination of preselection criteria can affect the 

Chapman et alll200l 



derived redshift distribution 



I Lindner et al.|[20Tlt lYounger et al 



120091120071) ."and the 



need for optical spectroscopy biases against lensed sys- 
tems for which the optical redshift will correspond to the 
foreground galaxy. Photometric redshifts obtained us- 
ing submm bands are very useful for estimating the high 
redshift nature of the submm sources, but suffer from er- 
rors due to the degeneracy between the dust temperature 
and the redshift , which limit their precision to Az ps 0.3 
(|Aretxaga et al.l 120071: IHughes et all 120021 ) . When the 
photometric redshift estimates involve SED template fit- 
ting, errors can also arise from our limited knowledge of 
the intrinsic SMG SED from FIR to radio, and its evo- 
lution with redshift. Direct spectroscopic redshift deter- 
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mination at submm wavelengths allows us to avoid such 
problems and study SMGs over a wide range of redshifts. 
Combined with multi-wavelength data, training sets of 
spectroscopic redshifts may also prove useful for reduc- 
ing these errors for application to the large photometric 
datascts from ongoing and future surveys. 

SMGs con tain large reservo irs of molecular gas 
(10 10 - 11 M , ITacconi et aTll2008D . whose cooling is dom- 
inated by the rotational lines of CO, almost equally 
spaced by ~115 GHz in the rest frame. Thus, the CO 
line detections at wavelengths beween 1 cm and 1 mm 
(30-300 GHz) offer the most direct measurement of their 
redshifts. Howeve r, with the exception of only three 



y three 
1 l2009at 



other CO redshifts (iDaddi et al.l l200a IWeifi et al 
iSwinbank et al.ll2010D . prior to the Herschel surveys the 
CO detections have largely been limited to SMGs whose 
redshi fts were already kn own from optical spectroscopy 
(e.g., iFraver et al.lll998T ). as a consequence of the small 
instantaneous bandwidth of typical mm-wavelength re- 
ceivers. This picture is rapidly changing with the ad- 
vent of a new gener ation of instrum ents, such as Z-Spe c 
(|Navlor et al.ll2003t ). Zpectrometer (jHarris et al.1 120071 ). 
and the new receivers used on interferometers such as the 
IRAM Plateau de Bure Interferometer (IRAM/PdBI), 
the Combined Array for Research in Millimeter- wave As- 
tronomy (CARMA), and the Atacama Large Millimeter 
Array (ALMA). Z-Spcc overcomes the mentioned lim- 
itations due to its large bandwidth, covering the en- 
tire 1-1.5 mm atmospheric window, which allows simul- 
taneous observations of multiple CO lines for galaxies 
at redshifts z > 0.5. Although the potential of using 
the C O ladder for redshift determination is we ll known 
(e.g., iCombes et aTlll999t [Sanders et al.lll986f ). due to 
sensitivity limitations of current instruments, only large 
area submm surveys can provide a significant number of 
sources bright enough for such measurements. 

These spectra can be used not only for an efficient 
redshift determination, but also to constrain the phys- 
ical properties of the gas and dust (e. g., mass, den- 
sity, t emperature) in these galaxies (e.g., [Bradfor d et al.1 
I2009T ) . by measuring the CO line strengths and the con- 
tinuum slope. The analysis of the CO properties re- 
quires measurements of multiple CO lines, often involv- 
ing the use of multiple instruments. To date, several 
spectral line energy distributions (SLEDs) for the CO 
molecule have been constru cted for small mixed sam- 
ples of galaxies and quasars (iPapadopoulos et al.ll2010bt 
iWang et al.l 120101 iBavet et all 120091 ) . or individual ob- 
jects. Relatively well sampled CO SLEDs have been 
constructed from the ground for some bright quasars 
(jWeifi et al.ll2007atlBradford et al.ll2009ft . while complete 
CO SLEDs have been measured by the Herschel Space 
Obse r vatory in low red s hift g alaxies (jPanuzzo et al.1 
l20T0l Ivan der Werf et al.l I27M . Most SMGs have 
been observ e d in o nly one or two CO l ines (see e.g., 
Harris et ail 120101: llvison et al.l l2010at ITacconi et al.l 
2008 HGreve et al.ll2005HSolomon fc Vanden Bouti^lM ? 
and their physical properties remain largely unknown. 
This situation has improved in recent years, with ob- 
servations of multiple CO lines in individual SMGs 
(lAo et al.1120081: iCarilli et al.ll2010l: ILestrade et al.l 120101: 



Riechers et al.l 20101 : iDanielson et al. 1 120101: IScott et al.l 
2011bl) . Thebest sampled CO SLEDs show that multiple 



CO components are required to explain the full line lumi- 



nosity distribution, where most of the mid- J CO emission 
can generally be fit by a warm component, with kinetic 
temperatures of 40-60 K and gas volume densities of 10 3 - 
10 cm -3 . However, solutions with kinetic temperatures 
of a fewx 100 K and lower densities are also allowed by 
the d ata (|Ao et al.ll20"ol IWeifi et al.ll2007at IBavet et al.l 
|2009| ). and this region of the parameter space has been 
insufficiently explored. With Z-Spec we can cover some 
portion of the CO SLED in a single observation (depend- 
ing on the redshift), with a common calibration for the 
entire bandpass, and we can start to place broad con- 
straints on the parameter space. However, additional CO 
line measurements, especially for the CO(l— 0) line, can 
prove essential in distinguishing between possible mod- 
els, or identifying a substantial amount of cold gas. 

This paper describes observations of five H-ATLAS 
sources undertaken with Z-Spec. Based on the CO emis- 
sion detected by Z-Spec, we successfully determined the 
redshifts of four out of five targets, helping confirm that 
they are lensed. The Z-Spec observations are described 
in Section [2] followed by the description of the algorithm 
for redshift determination in Section [3] We use the mea- 
sured redshift to constrain the spectral energy distribu- 
tion (SED) of these galaxies, estimating the dust temper- 
ature and emissivity index, as well as the total infrared 
luminosity. We perform an analysis of the partial CO 
SLEDs, constructed from the lines observed by Z-Spcc, 
to constrain the physical conditions of the molecular gas. 
The analysis of the galaxy SEDs and CO emission lines 
is presented in Section 01 and a summary of our results 
can be found in Section [5] Throughout the paper we as- 
sume a standard ACDM cos mology, with Hp =71 km s _1 
Mpc-\ O M =0.27, O\=0.73 (|Spergel et al.l[2007h . 

2. OBSERVATIONS AND DATA REDUCTION 

We selected five high-z candidates among submm- 
bright galaxies with F(500/^m)>100 mJy (TablcQ]) from 
the H-ATLAS survey for follow-up observations with Z- 
Spec on the 10-m Caltech Submillimcter Observatory 
(CSO). The flux limit was chose n based on theoretical 
calculations (jNegrello et al.ll2"007D . which show that high 
redshift galaxies may have observed 500 fim fluxes above 
the 100 mJy threshold only if lensed by foreground ob- 
jects. In the H-ATLAS SDP catalogue 11 objects satisfy- 
ing the flux cut have been found, out of which 6 objects 
have been identified as contaminants (four nearby spi- 
rals, one Galactic s tar forming region, and one blazar, 
INegrello et al.ll2~010D . resulting in a total of five remain- 
ing lens candidates. For convenience, throughout the 
paper we identify our targets by their names used in 
the SDP H-ATLAS catalogue (SDP.9, SDP.ll, SDP.17, 
SDP.81, and SDP.130). In order to distinguish these 
submm-bright lens candidates from the foreground lens- 
ing galaxies, it was necessary to measure their redshifts 
directly at submm wavelengths and confirm that they 
are at higher redshifts than the foreground galaxies. The 
redshifts of the foreground objects have been separately 
measured in the optica l and near-infrared, and found to 
be in the range 0.3-0.9 (jNegrello et al.|[2010l ). much lower 
than the redshifts of the submm galaxies, thus support- 
ing the lcnsing scenario. Several instruments were in- 
volved in the submm redshift-dctermination follow-up: 
CSO/Z-Spec, GBT/Zpectrometer, and IRAM/PdBI. 
The GBT/Zpectrometer results have been presented in 
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Table 1 

Summary of the Z-Spec observations on the H-ATLAS sources. 



IAU Name 


H-ATLAS 


Dates 


T225GHz 


Integration 


rms Uncertainty^' 




SDP ID 


Observed 


(zenith) 


Time (hrs) 


(mjy) 



H-ATLAS J090740. 0-004200 
H-ATLAS J 091043. 1-000322 
H-ATLAS J090302.9-014128 
H-ATLAS J090311. 6+003905 
H-ATLAS J091304.9-005344 



SDP.9 Apr 27 - May 14 0.05 

SDP.ll Apr 28 -May 4 0.05 

SDP.17 Mar 28 - Apr 1 0.04 

SDP.81 Mar 7 - Mar 12 0.02 

SDP.130 Mar 21 - Mar 22 0.04 



0.21 10.6 4.0 

0.18 6.8 5.5 

0.08 18.2 2.9 

0.05 22.5 2.3 

0.08 8.6 4.4 



Note. — The columns list: 1) the IAU source identification; 2) the ID of the source in the SDP H-ATLAS 
catalogue; 3) the range of dates for the observations; 4) the range in T225GHZ over all observations of the source; 
5) the total integration time on the source (including the time spent in the off-source position during the nod 
cycle, but excluding all other overheads); and 6) the median rms uncertainty on the measured flux density. 

( a ' Varies with frequency. The channel width is frequency dependent, with a mean of 950 km s~ 1 . 
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Figure 1. The Z-Spec spectra of four submillimetcr bright H-ATLAS galaxies. The fit to the continuum and CO lines at the measured 
rcdshift is ovcrplotted in red, and the positions of the strongest lines falling in the Z-Spec bandpass are indicated by the vertical blue lines. 
The line indicated in red in the spectrum of SDP.130 is unidentified. 



iFraver et all (|2011l) . while this paper shows the CSO/Z- 
Spec results. 

Z-Spec is a single spatial pixel grating spectrome- 
ter with 160 silicon-nitride micro-mesh bolometer de- 
tectors (i.e. chann e ls) operating from 190 — 308 GHz 
(iNavlor et all IzOOl lEarle et all [200l Bradford etaLI 
12009|). The frequency response of the Z-Spec channels 
is approximately gaussian, with a variable FWHM from 
720 to 1290 km s^ 1 over the ba ndpass, that is ro ughly 
equal to the channel separation (Earle et al. 2006). The 
Z-Spec beam size FWHM at the CSO has been measured 
to vary from 39 to 25 arc-sec across the band. 



We carried out the Z-Spec observations of H-ATLAS 
sources at the CSO from 2010 March 07 - May 14 un- 
der generally good to excellent observing conditions, ac- 
cumulating from 6.8 to 22.5 hour integrations on each 
target. The zenith opacity at 225 GHz (monitored by 
the CSO tau meter) was t 225GHz = 0.06 on average, and 
7"225GHz < 0.07 for 75% of the observations. A summary 
of the observations, including the total integration time 
on each source is given in Table [1] The Z-Spec data were 
taken using the standard "chop-and-nod" mode in order 
to estimate and subtract the atmospheric signal from the 
raw data. The secondary mirror was chopped on- and off- 
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Figure 2. The Z-Spcc spectrum of the H-ATLAS source SDP.17. The fit to the continuum and CO lines at z = 0.94 is ovcrplotted in 
red in the upper panel, and the rotational CO lines are indicated by the vertical blue lines. These lines have been subtracted from the 
spectrum shown in the lower panel. The red line in the lower panel shows the fit including the lines identified at z = 2.308. 



source at a rate of 1.6 Hz with a chop throw of 90 arc-sec 
while stepping through a 4-part nod cycle which posi- 
tion switches the primary mirror, integrating for 20 sec 
at each nod position. The chopping removes atmosphere 
fluctuations and the nodding removes instrumental off- 
sets due to imperfect match between the two chopped 
positions. We checked the pointing every 2 — 4 hrs by 
observing quasars and other bright targets located close 
in elevation to the H-ATLAS targets, making small (typ- 
ically < 10 arc-sec) adjustments to the telescope pointing 
model in real-time. 

We analyze the data using customized software i n the 
same manner as described in Bradford efal] (|2009fl . For 
each channel, the nods are calibrated and averaged to- 
gether, weighting by the inverse variance of the detec- 
tor noise. Absolute calibration is determined by obser- 
vations of Mars once per night, which we use to build 
a model of the flux conversion factor (from instrument 
Volts to Jy) as a fun ction of each detector 's mean operat- 
ing ("DC") voltage (IBradford et alS^M i. Since the DC 
voltage depends on the combination of the bath temper- 
ature and the total optical loading on the detectors, we 
use these curves to determine appropriate calibration fac- 
tors to apply to each nod individually. Based on the root 
mean square (rms) deviations of the Mars measurements 
from the best-fit curves, the channel calibration uncer- 
tainties are 3 — 8%, excluding the lowest frequencies for 
which a clean subtraction of the atmosphere is hindered 
by the pressure-broadened 183 GHz atmospheric water 
line. These uncertainties are propagated through the 
data reduction. The median rms uncertainties on the fi- 



nal co-added spectra for the H-ATLAS galaxies are listed 
in Table HJ These errors do not include the ~ 5 % uncer- 
tainty on the brightness temperature of Mars (jWrightJ 
120071 ). The calibrated Z-Spec spectra of the five ATLAS 
galaxies are shown in Figures [T] and O The rcdshifts of 
these sources arc determined using a custom algorithm, 
tailored specifically for multiple lines observed simultane- 
ously in the same bandpass. This algorithm is presented 
in the next section. 

3. REDSHIFT DETERMINATION 

3.1. Algorithm Description 

The redshift determination relies on multiple CO and 
atomic lines being present in the Z-Spec bandpass. Since 
not all these lines are necessarily strong enough to be 
individually detected at high significance, we developed 
a redshift-finding algorithm that is capable of handling 
cases where the signal-to-noise in individual lines is low, 
by combining the significance of all these lines. The num- 
ber of CO lines redshifted in the Z-Spec bandpass grows 
from 2atz = 0.51 (CO(3-2) and CO(4-3)) to 4 or 
more at z > 2 (starting at CO(5-4) through CO(8-7)). 
Since under most excitation conditions present in Ultra 
Luminous Infrared Galaxies (ULIRGs) and SMGs the in- 
tensity of the CO ladder drops beyond ~CO(7— 6), it can 
become increasingly difficult to measure redshifts higher 
than ~ 3.2 in the absence of high-excitation, warm CO 
gas. 

For the redshift determination we use a reference 
line list containing the lines expected to be strong in 
ULIRGs and SMGs, namely the CO rotational lines (up 
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to CO(17-16)), the [C I] 492.16 GHz line, the [N II] 
1458.8 GHz line, and the [C II] 1900.569 GHz line. As 
the width of the Z-Spcc channels varies from 720 to 
1290 km s^ 1 over the bandpass, larger than most ob- 
served line widths, most of the signal from one line will 
be concentrated in a single channel. Therefore, in order 
to determine which lines are present in the spectrum, we 
need to look at the signal-to-noise in individual channels. 
The [C I] 809.342 GHz and CO(7-6) 806.651 GHz lines 
are blended in the same channel and therefore degenerate 
for the purpose of this procedure. Recent Herschel ob- 
servations suggest that water li nes might also be bright i n 
certain ULIRGs, like Mrk 231 (|van der Werf et al.|[20Toh . 
while the con firmation of the wa ter line in SDP.17b by 
IRAM/PdBI (lOmont et all 120 111 this paper, ) indicates 
that this might also be the case for some SMGs (sec 
Section [3. 4j) . However, we defer using such lines in a sys- 
tematic way until more data is available on the presence 
of water emission in ULIRGs and at high redshift. 

The significance of the determined redshift is depen- 
dent on which lines are present in the spectrum relative 
to the lines that were expected to be observed, based 
on the reference list (defined above). This is the case 
for SDP.17b (see Section l3T4| . where the redshift signif- 
icance increases greatly if we include the water line on 
our line list. However, such an extension of the refer- 
ence line list is not always justified, since the significance 
of the redshift of another galaxy where the water line 
is not detected may be unnecessarily diminished. Care 
must also be taken in using the current line list at higher 
redshifts, where the high-J CO lines are likely to have 
much lower significance relative to the [N II] and [C II] 
lines. The least biased way to introduce this constraint 
may be to require that [C II] be the brightest line in the 
spectrum, and/or to limit the range of CO lines searched 
for to lower J's. 

No a priori knowledge of the relative line strengths 
is assumed, and therefore the algorithm gives equal 
weight to all the lines in the reference list. Even 
though it is known that the strength of the CO 
lines drops with increasin g J for starburst galaxies 
(e.g., iDanielson et al.l[2010h . but remains relatively con- 
stant for AGN-dominated galaxies and qua sars (e.g., 
Bradford et "al|[200l Ivan der Werf et alj|2010j) . it is gen- 
erally impossible to know apriori the nature of the emis- 
sion in the galaxy being observed. Associating line 
weights according to a model might artificially increase 
the significance of some redshifts and decrease the signifi- 
cance of others. Moreover, this would not prevent a non- 
detection in the cases where the signal-to-noise does not 
pass our threshold criterion (for example SDP.130, Sec- 
tion I33J. We can always use such relative line strength 
templates as a consistency check for the redshift determi- 
nation, rather than as an integral part of the algorithm. 

The redshift finding algorithm uses two test statistics, 
E\{z) and E 2 (z) (Eqs. [T]and[2]), constructed from com- 
binations of the detection significance in those channels 
in which a reference line would be observed by Z-Spcc 
at redshift z. The values of these test statistics are re- 
lated to the probability that the lines from the reference 
list, redshiftcd by a factor of (1 + z), are present in the 
spectrum. 

Let N(z) be the number of reference lines that would 
fall in the Z-Spec bandpass at redshift z. We search a 



redshift range between 0.5 and 6.0 in steps of 0.001. How- 
ever, the redshift determination and the false detection 
rate are not sensitive to the exact redshift range being 
searched, as long as multiple lines fall in the bandpass 
and the actual redshift is included in the search. The 
algorithm loops through all the z values, redshifting all 
the lines in the line list, and finding the set of N(z) Z- 
Spec channels corresponding to the lines in the bandpass 
for each individual redshift. The two test statistics, E\ 
and Ei , are evaluated for each redshift using the contin- 
uum subtracted signal Si and the noise <Ji in the set of 
N(z) channels determined in the previous step. The con- 
tinuum subtraction uses a fourth degree polynomial to 
better account for local smooth deviations from a power- 
law. 

The first test statistic, E\(z), is defined as the ratio of 
the total signal to the noise, summed only over the Z- 
Spec channels that correspond to a line in our list when 
redshifted to redshift z, 



E 1 (z) 



Si -Si 



(i) 



where the sum is taken from 1 to N(z), and Si and Oi 
are the signal and noise, respectively, for the channel 
corresponding to line i. 
The second test statistic, E 2 (z), is defined as 



E 2 (z) = median{/y | /jj = 0.5(Si/cn + Sj/aj), 



(2) 



l<i,3< N(z),i< j}x ^/W{z) 
where the set contains all possible pairs of lines in the Z 



Spec bandpass at the corresponding redshift, and \/N(z) 
is a normalization factor, such that the distribution of 
E2(z) for a noise spectrum approaches a standard normal 
(A/"(0, 1), see Appendix). 
An alternate definition of E\ would be 



1 



\ Si 



(3) 



It can be shown (see Appendix) that for any individual 
redshift this estimator has a higher significance than E\ 
(larger expected value), which would make it a better 
choice when taken independently from E 2 . However, our 
simulations show that we obtain a lower number of false 
positives when using E\ rather than E$ in combination 
with E 2 , since E\ and E 2 are less correlated than E$ and 
E 2 . The details are given in the Appendix. 

All three statistics defined above are maximized when 
the redshiftcd frequencies of the reference lines match 
the frequencies of the channels with the highest contin- 
uum subtracted significance. Their distributions are well 
reproduced by standard normals for all redshifts when 
there is no signal in the spectrum (consistent with noise), 
since in this case all Si/ at have a standard normal dis- 
tribution (see Figure [3] and the Appendix). 

We consider a redshift secured when at this redshift 
both Ei and E 2 reach their maxima, and the signal-to- 
noise combination is larger than a certain threshold (de- 
fined in terms of a new statistic E 2max (zo) > 2.12; see 
Section I3.2[) . Even though the maximum of any of the 
two statistics could be used for redshift determination, 
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Figure 3. Distributions of the two test statistics derived from 
blank sky spectra. The histograms for E\ and E2 are shown as the 
black and dahsed red histograms, respectively. Gaussian fits corre- 
sponding to the listed standard deviations are overplotted in black 
and red, respectively. In the noise simulations, as well as in the 
sky spectrum, the E\ and E2 distributions will be well described 
by standard normals, since all Si/ai have also a standard normal 
distribution. 

the use of two statistics instead of one, as well as a signal- 
to-noise cut-off, helps reduce the number of false redshifts 
that can be due to random noise fluctuations in the spec- 
trum. For redshifts < 0.5, and possibly > 6.7, the pres- 
ence of only one line in the spectrum does not allow an 
unambiguous redshift determination. Note that E± (or 
equivalently E3) would be a reasonable statistic for single 
line detections, but note that E2 is undefined unless mul- 
tiple lines are present in the Z-Spec bandpass at a given 
redshift. The conditions for a secure redshift determi- 
nation when multiple lines are present in the spectrum, 
and the significance associated with the derived redshift 
are further discussed in the next Section. 

3.2. Noise Simulations 

In order to determine the properties of our estimators 
and the criteria for a redshift to be secured, we need 
to run noise simulations based on the actual measured 
Z-Spec noise in each channel, and construct the distribu- 
tions of these estimators. In the end, this will allow us to 
establish the significance of our redshift determination. 

The noise per channel is obtained from the power spec- 
tral density (PSD) of the time series for each nod. In the 
Fourier transform of the time series, the signal will be 
contained at the chopper frequency, and the noise is es- 
timated by averaging the values of the PSD around the 
chopper frequency. Our final co-added spectra contain 
nods from multiple observations, weighted by the indi- 
vidual noise estimates. The final uncertainty associated 
with the co-added spectra is calculated from the noise 
in all the ind ividual nods, following the prescription of 
IZhanal pOOl ) for weighted means. The calibration error 
is not taken into account because it affects equally the 
signal and the noise, leaving the significance per channel 
and the values of the test statistics unchanged, which is 
one of the strengths of this method. 

To be able to simulate the estimator behavior in the 
absence of any signal, we need to start by choosing a noise 
distribution. For our simulations, we assume that the 




Significance level (in 0.05 bins) Normal theoretical quantiles 

Figure 4. Left: Distribution of p-values for the K-S test, compar- 
ing the noise distribution Gik (Equation \4§ for each channel i to 
the standard normal, for all 160 channels. For p-values above 0.05 
we cannot reject the null at the 5% significance level. The dotted 
line shows the distribution of the K-S test probabilities for 160 sets 
of random points drawn from a A^(0, 1) distribution. Right: Q-Q 
plot of noise quantiles vs. standard normal quantiles, for a set of 
16 channels, color-coded from the lowest (indigo) to the highest 
(red) frequency. The corresponding K-S p-values are also color- 
coded. If the two distributions are identical, the points should lie 
along the diagonal. Overplotted with a dotted line is a simulated 
observation, with points randomly drawn from a A^(0, 1) distribu- 
tion, showing a scatter similar to our channels. The corresponding 
p-value is shown in black in the top left corner. 

noise is gaussian distributed, with a standard deviation 
given by the measured error in each channel. To test 
this assumption, we study the noise in a 6.5 hour Z-Spec 
integration on blank sky, recorded on May 11 and 12, 
2010. For each channel i for our blank-sky data, we look 
at the quantity 

G ik = (4) 

0~ik 

where fc represents the nod number, i is the channel num- 
ber, Ai is the average signal in channel i over all nods, 
and Sik and Oik are the signal and the noise, respectively, 
for the corresponding nod and channel. In our blank-sky 
data there are about 220 nods per channel, after flag- 
ging. The distribution of Gik for a given channel i, over 
all nods, should be a standard normal if our assumptions 
are correct. 

We first apply a Kolmogoro v-Smirnov test (K-S, 
iKolmogorovl |1933t iSmirnovl fl94l) . which tests the hy- 
pothesis that the observed noise distribution is drawn 
from a standard normal by comparing their cumulative 
distributions. The test results arc quantified in terms 
of the p-value, which is the probability that a value of 
the test statistic equal or greater than the one observed 
would be obtained if the null hypothesis were true. In 
the left panel of Figure HI we show a histogram of the 
p-values of the K-S statistic for all the channels, which 
demonstrates that we cannot reject the null hypothesis 
that the noise is gaussian distributed at the 5% level for 
any of the channels, since all p-values lie above this level. 
The plot also shows a large spread in p-values from chan- 
nel to channel. For comparison, we run the same K-S test 
for 160 sets of random numbers drawn from a 7V"(0, 1) dis- 
tribution. Each set contains the same number of samples 
as the corresponding channel. The distribution of the K- 
S p-values for these computer-generated normal samples 
is shown by the dotted line in the left panel of Figure |4j 
This simulation also shows a large spread in p-values, 
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approximately uniform across the range. 

The K-S statistic can be affected by a variety of factors, 
pertaining both to departures from gaussianity (shape of 
the distribution), and to mismatches in the parameters 
of the assumed distribution (i.e., the sample is drawn 
from Af(/xo,Oo) instead of Af(0, 1)). If our assumption 
that the noise for each channel is gaussian distributed is 
correct, but we have over- or under-estimated aik, this 
can, in principle, result in small p-values for the K-S 
test. A use ful tool in this case is the q uantile-quantilc 
(Q-Q) plot (|Gnanadesikan fc WilklH96l . which is more 
sensitive to multiple aspects of the distributions being 
compared, but does not provide a quantitative measure 
of these deviations. 

A Q-Q plot is basically a representation of the observed 
data quantiles versus the theoretical quantilcs of the as- 
sumed distribution (AT(0, 1) in this case). The quantilcs 
are defined as regular intervals on the cumulative distri- 
bution function, intuitively intervals of equal probability. 
The Q-Q plot is demonstrated for a sample of 16 channels 
in the right panel of Figure |4] If the noise is gaussian dis- 
tributed and the aik are estimated correctly, the points 
on Figure 0] for each channel will follow the diagonal. If 
the noise is over- or under-estimated, the relationship will 
be still linear, but with a different slope. Figure 0] shows 
that this might be the case for some of the low-frequency 
channels (blue) , for which the noise is known to be more 
variable due to both intrinsic bolometer problems and 
atmospheric noise. These very-low frequency channels 
are in fact excluded from our redshift-finding algorithm. 

Departures from gaussianity, which could be due to 
noise correlations between channels, would stand out in 
the Q-Q plot as departures from linearity. For compar- 
ison, we overplot in Figure 2] the curve obtained for a 
computer-generated sample drawn from JV(0, 1) (black 
dotted line) , which shows a similar level of scatter around 
the diagonal as our channels. We conclude that the scat- 
ter of the noise distribution around the linear correlation 
in the Q-Q plot is negligible when compared to the re- 
sults from the random number generator for a standard 
normal, supporting the results of the K-S test for the 
gaussianity of the channel noise. 

Our simulations create multiple realizations of a pure- 
noise spectrum, with the signal in each channel being 
a gaussian random variable with mean and standard 
deviation equal to the measured noise in that channel, 
A/"(0,af). Even if for some channels the noise might be 
over- or under-estimated, its exact value is not essential, 
since it cancels out as part of the signal-to-noise ratios 
in estimator definitions, and in the end we are left with 
A/"(0, 1) distributions for the estimators (see Appendix). 
As discussed above, the noise per channel from our blank- 
sky observation is well approximated by a gaussian distri- 
bution, aside from small channcl-to-channel noise corre- 
lations. In our simulations, we reproduce these residual 
noise correlations between different channels using the 
method of Cholcsky factorization. The noise correlation 
matrix is constructed from all the nods contained in the 
blank sky spectrum, 

r EkiSik A)(S jk - Aj) 

v sqrfZ k {S ik -A i YY, k {S i k-A j r U 
where the sums are taken over all nods. After multiply- 



ing a randomly-generated uncorrelated vector with the 
lower-triangular matrix from the Cholcsky decomposi- 
tion, one obtains a vector with the sa me correlation prop- 
erties as the original sky noise model (Kaiser fc DickmarJ 

Eia e.g., )■ 

We run separate simulations for correlated and un- 
correlated noise, each with 10 5 realizations of noise spec- 
tra. For each realization of the noise spectrum, we record 
the maximum values of the estimators E\ and E2, and 
construct their joint distribution function over all real- 
izations. For any measured max(Ei) and max[Ei) from 
real data, we define the associated false detection rate 
(FDR) as the probability of finding a maximum value 
of Ei > max (Ei) and Ei > max(E^i) by chance, in the 
absence of real signal. We calculate this joint probabil- 
ity from the simulated 2-D right-cumulative distribution 
function of the maxima of the two estimators, as shown 
in Figure[7]for each measured (max(Ei), moif^)) pair. 
Figure [H] shows the marginal FDR (the 2-D cumulative 
distribution marginalized over E2) for all the max(Ei) 
values (solid black line), as well as for the max(Ei) val- 
ues left after imposing the additional constraints on the 
estimators discussed below. 

As the first constraint, for each simulated spectrum we 
identify all the max(Ei) and max(E2) values that satisfy 
the condition that both estimators reach their maxima at 
the same rcdshift. As can be seen from the dashed line in 
Figure [51 imposing this condition reduces the total FDR 
to about 40% for both correlated and un-correlated noise. 
The decrease in FDR is due to the fact that the locations 
of the maxima of Ei and E2 deviate from a perfect corre- 
lation (Pearson correlation coefficient < 0.80), with the 
scatter spread over all redshifts. This shows that the 
combination of two test statistics is more robust against 
random fluctuations than any estimator used indepen- 
dently, and considerably reduces the noise floor across 
the redshift range. 

Requiring that the two estimators be maximized at the 
same redshift we still get a rather high total FDR (at 
least 40%). In order to further reduce the number of 
spurious redshifts obtained from blank sky spectra, we 
introduce a signal-to-noise threshold cut. We define the 
quantity 

E2max(z) = max{/y \fij = 0.5(S t /a t + Sj/aj), 
l<i,j<N(z),i<j}. 

This definition is very similar to E2(z), with the median 
replaced by the max, and without the normalization fac- 
tor. The normalization factor is not needed here be- 
cause we want to be able to establish a threshold crite- 
rion across the entire redshift range, independent of the 
number of lines N(z). Since this is not an estimator, 
we are not interested in standardizing its distribution, 
and moreover, its distribution will be likely different from 
that of E2, as an extreme order statistic rather than a 
central order statistic. 

For any i and j the quantity /y = 0.5(Si/o~i + Sj/aj) 
is distributed as a normal with standard deviation o~ av = 
l/y/2 (since Si/o~i has a standard normal distribution). 
We choose the value £2 mM (z) = 3a av = 2.12, as our 
threshold. In other words, we require that strongest pair 
of lines at the determined rcdshift have an average signal- 
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Figure 5. Null test for the redshift finding algorithm, using a 
blank sky spectrum. The dashed horizontal lines are drawn at 
and 3, to guide the eye. Note that the values of the estimators are 
below 3 across the redshift range, the positions of their maxima do 
not coincide, and the E^max values corresponding to these maxima 
are below our threshold. The dashed vertical lines indicate the 
positions z(max(Ei), and z(max(E2)), respectively. 

to-noise greater than 2.12. This places approximate lim- 
its on the signal-to-noise per channel and the values of 
the estimators for the determined redshift of ~2 and ~3, 
respectively. If no signal is present in the spectrum (null 
hypothesis), the distributions of E\{z) and E 2 (z) are well 
approximated by standard normals for all redshifts z (see 
Figure [3]) . As such, our noise cut also implies rejecting 
the null hypothesis at the ^99% level for the determined 
redshift. A realization of the two test statistics using the 
blank sky spectrum is shown in Figure[5J In this case, the 
maxima of the two statistics occur for different redshifts, 
the values of E^max are below our threshold for both E\ 
and E2 maxima, and the values of both estimators are 
below 3 for all redshifts. 

The effect of imposing the Ei max constraint on the 
FDR is shown by the dot-dashed curve in Figure [51 
The noise threshold criterion improves the significance 
(1-FDR) of the lowest signal-to-noise results by cutting 
the limiting FDR down to about 50% for un-correlated 
noise and to ~ 75% for correlated noise. The joint effect 
of both constraints is shown by the dotted curve in Fig- 
ure |H1 and the corresponding 2-D FDR distributions as 
a function of the maxima of both E\ and E2 are plotted 
in Figure [7] We derive that the estimator values pass- 
ing these 2 tests will result in a total false detection rate 
lower than 24% for uncorrelated noise and below 33% for 
correlated noise. 

We emphasize that these FDR values are limiting val- 
ues, in the sense that they are independent of the actual 
values of max(Ei) and max{E2), and only satisfy the 
requirement that the estimators pass the two tests (at- 
tain maxima at the same redshift and exceed the E2max 
threshold). For any actual redshift determination, the 
associated false detection rate will be determined by the 
values of max(Ei) and max(E2) for that particular spec- 
trum, which are inversely correlated with the FDR, as 
indicated by the vertical lines in Figure El For the pri- 
mary redshifts determined in this paper, we find that 
the integrated false detection rates are smaller than 2%, 
as described in Section 13.41 and listed in Table O Note 
that all quoted FDRs arc by definition integrated over 
the whole redshift range, and represent the probability 
of obtaining a false positive, at any redshift, given the 



values of the estimator maxima, but the probability of 
obtaining the same redshift by chance is roughly a fac- 
tor of ~5000 lower, based on the number of redshift bins 
searched. This is because it is even more unlikely that 
the same channels fluctuate high due to noise alone. 

To summarize, a redshift zq is accepted when the fol- 
lowing conditions hold: 

Ex(zq) = max(E'i), 

E 2 (z ) = max(E 2 ), (7) 

E2max(zo) > 2.12, 

where the last condition is basically a signal-to-noise 
threshold criterion, and the significance of the estimated 
redshift is calculated as 1-FDR, with the FDR derived 
from the noise simulations, as explained above. 

At the CSO, Z-Spec can reach a measured maximum 
sensitivity of 0.5 Jy s 1 / 2 p er channel for an a tmospheric 
optical depth t 22 5=0.068 (jlnami et al.ll2008| ). Combin- 
ing the signal-to-noise threshold criterion with the mea- 
sured sensitivity of Z-Spec, we estimate that a redshift 
can be determined in less than 1.4 hours of integration 
time if the line flux densities per channel are on the or- 
der of 15 mjy, but can require more than 12.6 hours 
if the flux density is less than 5 mjy. For our galaxy 
sample, the mean integrated CO line flux (Table |3|) is 
~18 Jy km s _1 , while the average width of the channels 
is 950 km s . However, the flux density per channel 
could be only ~10 mjy if the line flux happens to be 
split between two adjacent channels. In this case, the 
typical integration time for obtaining a redshift with Z- 
Spec would be at least 3.5 hours. These time estimates 
reflect closely the best performance of the instrument 
and do not include calibration overheads. The actual in- 
tegration time needed to obtain a redshift will depend 
strongly on the instrument sensitivity at the time of the 
observations. 

3.3. Redshifts for the H-ATLAS SDP Sample 

The results of applying this algorithm to our galaxy 
sample is shown in Figures [8] and [7l We secure the red- 
shifts for four out of five sources, with an FDR < 10% 
in all cases. The redshift value and its uncertainty (Ta- 
ble [2]) are determined from the position and width of the 
peak of the E\ test statistic (Figures [8]). Due to the fi- 
nite width of the spectral channels, nearby redshifts can 
have the same or similar significance, since the lines will 
fall on the same channels for a narrow range of redshifts, 
given by our redshift space sampling. As we go further 
from the real redshift, some of the lines might still fall on 
the same channels, but not all of them, so the value of 
Ei will drop. We fit a gaussian to the peak of E\(z), and 
define the redshift error as the upper limit for the stan- 
dard deviation of this gaussian. This value is at least as 
large as the channel width, and accounts for the varying 
channel widths across the bandpass. 

For all the galaxies except SDP.130, the maxima of 
E\ and E2 satisfy both our criteria for a secure redshift 
determination, and for SDP. 17 we can identify a second 
redshift satisfying our criteria after subtracting the first 
set of lines from the spectrum. We calculate the signif- 
icance of the redshift for each of our sources by inter- 
polating the FDR at the observed values of max(E 1 ) 
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Table 2 

Summary of the H-ATLAS galaxy sample and the parameters derived from fitting their submm SEDs. 



H-ATLAS 


/< 


z 


Significance 






a 






At SFR 


SDP ID 






(%) 


(10 13 L ) 


(K) 




(10 9 M ) 


(arcsec 2 ) 


(10 3 M /yr) 


SDP.9 




1.577±0.008 


100 (99.97) 


4.4±0.5 


57±1 


3.8±0.2 


2.5 


0.65 


6.6±0.8 


SDP. 11 




1.786±0.005 


99.98 (99.22) 


7.8±0.9 


69±1 


5.7±0.4 


1.7 


0.43 


11.7±1.3 


SDP.17a( a > 




0.942±0.004 


87.33 (74.77) 


0.4±0.09 


27±1 


2.9±0.1 


4.9 


1.44 


0.6±0.1 


SDP.17b( a > 




2.308±0.011 


99.86 (97.46) 


3.9±0.9 


66±1 


2.9±0.1 


1.1 


0.30 


5.8±1.5 


SDP.81 


18-31< b ) 


3.037±0.010 


98.26 (90.02) 


6.4±0.3 


58±1 


3.2±0.1 


2.2 


0.69 


9.6±0.4 


SDP. 130 


5-7< b ) 


2.626±0.0003< c > 


N/A 


4.3±0.2 


55±1 


2.7±0.3 


1.6 


0.47 


6.5±0.3 



Note. — The columns list: 1) the ID of the source in the SDP H-ATLAS catalogue; 2) the gravitational lensing magnifica tion 
factor; 3) the measured redshift; 4) the redshift significance, calculated as 1 — FDR, where the FDR has been defined in Section 13.21 
The significance for correlated noise is given in parenthesis; 5) the integrated IR luminosity, obtained as the average between the 
SED fits with the CE01 libraries and the DH02 libraries. The factor fi is shown in front of quantities affected by gravitational lensing 
magnification; 6) the dust temperature, with the caveats described in the text; 7) the index of the power-law continuum fit to Z-Spcc 
data; 8) the dust mass; 9) the solid angle subtended by the dust emitting region; and 10) the star formation rate. 

( a ) The total observed flux was split between the two components, using a frequency-independent scale factor. 

( b ) Values taken from INegrello et al.l l20TOh . 

( c ) Redshift determined by GBT/Zpectrometer, followed by a more precise measurement with PdBI/IRAM (Ncgrcllo ct al. 2010). 

( d ) The uncertainties for Td are likely underestimated. The values shown are formal errors from the fit, and don't include correlations 
between parameters, or the inaccuracy of the assumed shape of the SED model. The values for /3 and vq arc kept fixed for all sources. 
( e ' Calculated in the optically thin limit. The dust masses should be interpreted as robust lower limits for the true total dust mass 
in the galaxy (see text). 
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Figure 6. FDR distributions resulting from noise simulations with un-correlated and correlated noise (left and right, respectively). The 
curves show the decrease in FDR as a function of the E\ value, before (solid) and after applying the condition E2 ma x > 2.12 (dash-dot), 
z(max(E\)) = 2 (max (£2)) (dashed), and both conditions jointly (dotted). The solid vertical lines correspond to the E\ values at the 
determined redshift for each target. These redshifts satisfy both conditions (dotted line), and the corresponding FDR for each of them is 
better shown in Figurc[7] For SDP.130 the values of the estimators at z=2.626 (GBT/Zpectrometer) are negative, and therefore cannot be 
maxima (no corresponding vertical line). The stars show the FDR associated with the values of max(Ei) obtained after subtracting the 
lines from each spectrum (the spectrum of SDP.130 is used as-is, and for SDP.17 after subtracting both sets of lines), each of them being 
placed on the curve corresponding to the conditions satisfied by the residual maximum. After subtracting the lines, none of the maxima 
satisfies all our criteria for redshift determination, and the individual significance for each of them can be read off these curves. 



and max{E-2). Figures [7] and |B] show the derived estima- 
tor values for each source relative to the FDR distribu- 
tion, and the significance of each redshift determination 
can be read off directly from these figures. Of all red- 
shifts that passed our criteria, the redshift of SDP. 17a 
(the second derived for SDP.17) has the lowest signifi- 
cance, close to 90% (75% for correlated noise), while all 
primary redshifts have a significance of at least ~ 99% 
(90% for correlated noise), equivalent to ~ 3cr or greater 
for a gaussian distribution. 

When lines are present in the spectrum, aside from 
the main peak due to the true redshift, secondary peaks 



will arise in the E\ and E 2 distributions, correspond- 
ing to redshifts where some of the lines in the line list 
fall on the same channels as the observed lines. The 
secondary peaks arc marked by blue asterisks for each 
source in Figure [5J The real redshift will have higher sig- 
nificance than the redshifts corresponding to these sec- 
ondary peaks, since the largest number of lines add their 
contribution to the total signal in this case. 

After removing the lines corresponding to the mea- 
sured redshifts from the spectra (both redshifts for 
SDP.17), the secondary peaks in the E\ and E2 distribu- 
tions are reduced to the noise level, and the maxima of 
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Figure 7. FDR contour plots from our simulations, as a function of E\ and E2 values. The points corresponding to the maxima of E\ 
and E2 for each of the sources that pass the cuts (all except SDP.130, see Figurc[8]l are shown by the filled stars. The contours correspond 
to the dotted line in Figure [6] (satisfying both z(max(Ei)) = z(max(E2)) and E^max > 2.12 conditions). The points obtained after line 
subtraction do not satisfy either condition and therefore are not shown. 



Ei and E2 fail to satisfy one or or both of our criteria. 
The FDRs associated with these line-subtracted spectra 
are plotted as stars in Figure [SJ indicating the signifi- 
cance of the remaining features. The stars are vertically 
positioned on the false detection curve corresponding to 
the criteria satisfied by max(Ei) after removing the lines. 
The computed FDRs and the fact that they do not sat- 
isfy both criteria indicate that in all cases the results 
after line subtraction are consistent with noise. 

All the measured redshifts, with their error bars and 
associated significance (calculated as 1 — FDR) for both 
correlated and un-correlated noise, are listed in Table El 
The significance of our redshift determinations together 
with the statistical redshift determination criteria (Eq.[7]) 
show that in general a redshift was already secured by 
Z-Spec after an integration time much shorter than the 
total integration time listed in Table [T] Future submm 
instruments with better sensitivity will be able to obtain 
the redshifts of such galaxies even faster, and open the 
possibility of large submm redshift surveys. 

3.4. Comments on Individual Redshifts 

The individual redshifts are presented in the order of 
the observations (see Table [1]). 

SDP.81 The redshift for SDP.81, z=3.037±0.01, ob- 
tained by this method on 19 March 2010, was confirmed 
(z=3.042±0.001) with follow-up observations with the 
IRAM Plateau de Bure Interferometer on 23 March 2010 
(IRAM/PdBI, INeerello et alJ20ialNerTeFall20Toh and 
with an independent blind search on 25 March 2010 
by the Zpectrometer instrument at the Green Bank 
Telescope (GBT/ Zpectrometer, INeerello et ail l2"oiot 
iFraver et al1l2011[ ). Both follow -ups were informed by a 
concurrent photometric redshift estimate (2.9+ ; 3 ). With 
the possible exception of the second redshift for SDP.17, 
this is the redshift with the lowest significance in our 
sample, with an E\ peak of 3.8, due to the weakness of 
the CO lines beyond CO(7 - 6). This is the first blind 
redshift obtained by Z-Spec. 

SDP.130 SDP.130 has a redshift of 2.6260±0.0003, 
measured by GBT/Zpectrometer (z=2.625±0.001, 



IFraver et al.1 20111). and made more precise with 
PdBI/IRAM (|Negrello et aTJIMol: iNeri et al.ll2010D . So 

far, three CO lines have been measured in this galaxy 
at this redshift, namely the CO(l— 0) line observed with 
the Zpectrometer, an d the CO(3— 2) and C O(5— 4) lines 
observed with PdBI (jNeerello et al.l 120101 ) . on a tuning 
that was successfully guided by the sub-mm photometric 
redshift of z = 2.6±g;| (INeerello et al J 12010ft . However, 
we do not detect any of the higher J transitions ( J u >6) 
that would fall in the Z-Spec bandpass at this redshift. 
The values of our estimators for z=2.625 are negative, 
suggesting that there is no signal left in the spectrum at 
this redshift after continuum subtraction. The estima- 
tors do not pass our redshift determination criteria for 
any other redshift, and Figure [6] shows the significance of 
the maximum E\ value obtained under these conditions 
(orange star). This non-detection, which places upper 
limits on the integrated fluxes of the CO(6— 5) through 
(9-8) lines of < 12.5 Jy km s _1 , suggests a low (< 50 K) 
gas temperature in the z=2.626 galaxy. We attempted to 
identify the line at 277 GHz, marked in red in Figure [1] 
with the CO(3— 2) transition at z=0.25, but that would 
be inconsistent with the optical sp ectroscopic redshift o f 
the lensing galaxy (0.220±0.002, INeerello et all 12010ft 
by more than 7000 km s _1 , as well as inconsistent with 
the observed SED. Based on the correlation beween 
the CO and the total infrared luminosity, the observed 
luminosity of the CO(3— 2) line would correspond to an 
ULIRG-class object at z=0.25, which would dominate 
the SED at 250 fim. No separate 250 /j,m-bright object 
is found nearby, and the PACS and SPIRE photometry 
of SDP.130 (see also Section 14. ip is inconsistent with 
the two sources being blended. Similarly, identifying 
this feature with the 987 GHz water line at z=2.626 
would require a velocity offset of ^4200 km s" 1 , and 
usually the presence of highly excited CO gas, which is 
not observed. This feature remains unidentified. 

SDP.11 Given the size of the Z-Spec beam 
(FWHM«30") and the possible presence of lensing or 
other foreground structures in the same beam, the ob- 
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Figure 8. Results of running the redshift-finding algorithm for all the H-ATLAS sources in our sample. The E2 test statistic has been 
offset vertically by 8 units, for clarity. The blue asterisks show the positions of the largest secondary peaks arising from coincidences with 
the lines from the actual redshift (see text). These peaks contain the same information as the main peak. In the SDP.17 panel, we note the 
extra peaks that do not match the secondary peaks corresponding to the first selected redshift. The SDP.17a panel shows the determination 
of the second redshift from the same spectrum, after subtracting the high-rcdshift component. No redshift is determined for SDP.130. 



served spectrum could be a combination of features from 
multiple objects. We choose this interpretation for the 
spectrum of SDP.17, best described by two components 
at different redshifts (both listed in Table [5]) . The first 
redshift found by our algorithm is 2.308 (SDP.17b). After 
fitting the CO lines at this redshift and subtracting them 
from the spectrum, we perform a second redshift determi- 
nation, identifying a second component with a redshift 
of 0.942 (SDP.17a). This combination explains all the 
features present in the spectrum (see Figure [2]), and is 
consistent with the interpretation of the 299 GHz feature 
as the restframe 987 GHz water line at a redshift of 2.308. 
This water line has been seen to be very strong in other 
AGN and star- formi ng galaxies at low redshift, su ch as 
Mrk231 and Arp 220 ([Gonzalez- Alfonso et all2010t) . and 
it has been t entatively detec t ed in the Cloverleaf quasar 
at z=2.56 bv Bradford et al.l (|2009D . More recently, mul- 
tiple excited water trasitions have also b een detected in 
the quasar APM 08279+5255 at z=3.91 dBradford et al.l 
l20TllLis et ^[2Wllvar7 der Werf et al.1120111) . The red- 



shifts of SDP.9 and SDP.17b have recently been con- 
firmed by follow-up observations of the CO(2— 1) and 
(3—2) lines, respectively (L. Leeuw, private communica- 
tion), with CARMA . Recent IRAM/PdBI observations 
(jOmont et al.1 120111 ) have confirmed the water line at 
z=2.3052, but did not find any other high significance 
line in the bandpass, which does not exclude the possi- 
bility of a CO(5-4) line at z> 0.944. The second red- 
shift (SDP.17a) has a much lower significance, but it is in 
agreement with the photometric and spectroscopic opti- 
cal redshifts (0 77±0 .13 and 0.9435±0.0009, respectively 
iNegrello et alJ[2~010( ). Alternatively, the peak now iden- 
tified with the CO(5— 4) line at z=0.94 could be arising 
from correlated noise fluctuations with the nearby water 
line. To confirm the presence of CO at z=0.94, we are 
planning a follow-up of the CO(4— 3) line. The presence 
of multiple ULIRGs in a single line of sight is intriguing, 
and is an example of discoveries that can be made possi- 
ble by Z-Spec's broad bandwidth. It also raises the pos- 
sibility that the flux-limited sample is affected by chance 
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Table 3 

Integrated fluxes for the emission lines identified in each galaxy. 



Line 




Frequency 
(GHz) 


SDR 9 


SDP.ll 


Integrated Line Flux (Jy km s ) 
SDP.17a 


SDP.17b 


SDP.81 


CO (4- 


■3) 


461.041 






14±9 






CO (5- 


-4 


576.268 


25±5 


23±8 


27±9( b ' 






CO (6- 


-5) 


691.473 


33±7 


29±10 




17±5 




CO (7- 


•6 


806.652 




18±14( a ' 




H±7(a) 


12±4( a ) 


CO (8- 


-7 


921.800 








16±6 


5±3 


CO (9- 


-8) 


1036.91 










6±3 


CO (10- 


-9) 


1151.99 










< 6.5 


[C I] 3 Pi - 


-> 3 P 


492.160 


< 11 




< 6 






[C I] 3 P> - 


^ 3 Pi 


809.342 




3 1±1 4(a) 




13±7 (a) 


< 6.5^ 


H a O 2 ,2 - 




987.914 








19±7(b) 





Note. — The columns list: 1) the transition; 2) the rest frame frequency of the transition; and 3) the inte- 
grated line flux for each galaxy (as measured, uncorrected for dust absorption) with 68% confidence intervals. 
Upper limits are 3<r. 

( a ) These lines originate in the same source and are blended at the Z-Spec resolution. The error bars account 
for this uncertainty. 

< b ) In the spectrum of SDP.17, the water line at z=2.308 and the CO(5-4) line at z=0.94 are blended. 



alignments, and the presence of multiple sources in the 
beam. However, this is likely a negligible effect for lensed 
sources, as the continuum sub-mm and mm flux will be 
clearly dominated by t he lensed, high-z gala xy, and not 
by the foreground lens (jNegrello et al.|[2010t) . 

SDP.9 and SDP.ll The significance of the redshifts 
for these galaxies corresponds to max(Ei) values of 6.5 
and 5.3, respectively (Figures [6] and fig:sig). The rcdshift 
of SDP.9 has been confirmed by CARMA observations, 
and more follow-up observations are currently planned 
for both SDP.9 and SDP.ll. 

4. GAS AND DUST PROPERTIES 

A model including the lines and power-law continuum 
is fit to each spectrum in Figures [T] and [3J allowing the 
line intensities, redshift, and continuum slope to vary. 
The best fit power-law index a for each galaxy is listed 
in Table [H The initial estimate for the redshift is pro- 
vided by the algorithm described above, and the fit is 
constrained by the requirement that all the lines be at 
the same redshift. In cases where some of the lines 
are blended, we first fit only the unblended lines to ob- 
tain a more precise value for the redshift, and then we 
fit all the lines simultaneously, with the redshift kept 
fixed, to get the integrated line strengths, listed in Ta- 
ble [3] Although the lines are not resolved, the signal 
from one line can be spread among adjacent channels 
due to the overlap of their frequency responses. We 
measure only the integrated line strengths, taking into 
account the frequency response of each Z-Spec chan- 
nel, weighted according to the line width. On average, 
line widths below ~1000 km s _1 are not resolved by 
Z-Spec, and we choose a value of 300 km s" 1 in fit- 
ting the integrated line strengths. This value closely 
matches the width of the lines for SDP.81 and SDP.130 
at PdBI (|Neri et al.ll2010f ). but is relatively low compared 
to the range found by interfcrometri c measurements of 
other lensed high-re dshift galaxies dGreve "eTaLl [20051: 
iKnudsen et all [2009] ). The CO(1-0) line widths deter- 
mined by the GBT/Zpcctromcter are somewhat larger 
( 435±54 km s" 1 for SDP.81 and 377±62 km s" 1 for 



SDP.130), suggesting that an additional gas component 
might contribute to this line. However, the determi- 
nation of the integrated line fluxes is not sensitive to 
the choice of the line width up to values of the order 
of the channel width. The largest uncertainties in the 
integrated line strengths arise in the case of line blend- 
ing, such as the CO(7-6) and [C I] 3 P 2 -^ 3 Pi lines, 
or the overlapping lines at different redshifts in SDP.17 
(blended lines are indicated in Table [3]). 

4.1. Continuum Spectral Energy Distributions 

The continuum data for all 5 galaxies is shown in 
Figure [9] The measured continuum flux from Z-Spcc 
is found to be in go od agreement with t he MAMBO 
1.2 mm photometry (jNegrello et al.l l2010( ). except for 
SDP.9. Estimates of the total amount of dust and star 
formation rates in each galaxy can be obtained by fit- 
ting their far-infrared (far-IR) to submillimcter spectral 
energy distribution (SED). In this fit we include the Z- 
Spcc data along with the Hcrschcl-SPIRE and Hcrschcl- 
PACS photometric points, as well as the Sub-millimeter 
Array (S MA) measurements a t 880 /xm for SDP.81 and 
SDP.130 (INeerello et al.l[20Tol) . 

The far-IR rest frame SED can be described by a mod- 
ified blackbody function, defined as 

F v = Q„(P)B v (T d )CL d 

- (1 r -TM("M f ') 2hvS 1 n, 

~ [ ' c ? e hv l kT * - 1 (8) 

L IR Q v {p)B v {T d ) 

4ird 2 jQ„(j3)B v (T d )dv' 

where Q v = 1 - e^" ^ '/ V °) P is the emissivity, B v (T d ) 
is the Planck function, r(z^o) = 1 is the optical depth 
at i/rj, fl d represents the observed solid angle of the dust 
emitting region, d is the (known) distance to the source, 
and h and k denote the Planck and Boltzmann constants, 
respectively. The fit can be performed with three param- 
eters: T d , /3, and a scale factor, while keeping vq constant. 
Including Vq as a fourth parameter in the fit leads to a 
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Figure 9. The best-fit SED models for the five H-ATLAS galaxies in our sample. The continuous line shows the modified blackbody 
spectrum with uq = 1300 GHz and /3 = 2.0, while the dotted and dashed lines show the SEDs obtained from the SED libraries of CE01 
and DH02, respectively. The total infrared luminosities are calculated as the average between the CE01 and DH02 SED template fits, to 
account for emission above the blackbody spectrum at higher frequencies. The parameters for the modified blackbody fits are also listed 
in Table [2] 



value of 1251±130 GHz for SDP.9, but no strong con- 
straints are found for the rest of the sample, leading us 
to fix the Vq at 1300 GHz. The low value found for Vq, 
and the observed flattening of the peak of the SEDs in 
the far-IR suggest that the SEDs of the galaxies in our 
sample can be modeled either as combinations of multi- 
ple graybodies with different temperatures, or as a sin- 
gle gray body with a large optical de pth at far-IR wave- 
lengths dPapadopoulos et al.ll2010H ). The overall scale 
of the SED can be parametrized either in terms of the 
solid angle fi^ or the total infrared luminosity (L/^), 
defined as the integral of the SED from 8 to 1000 /im 
(rest frame) . The Lj r derived in this manner underesti- 
mates the true total infrared luminosity, due to the likely 
presence of warmer dust components that contribute at 
shorter wavelengths. 

The simplest model that can r eproduce the data for 
the e ntire sample has fixed /3 = 2 (jPriddev &: McMahonl 
|2001[ ) and the already-mentioned ^o=1300 GHz, in agree- 
ment to the value found for SDP.9. The dust emissivity 
index /3 = 2 is also consistent with the error bars of the Z- 



Spec spectra. The best- fit models are shown in Figure GO 
and the corresponding values for Td are listed in Table [5J 
With dust temperatures between 54 and 69 K, the peak 
of the restframe dust SED is found in a narrow range of 
wavelengths (73 to 92 /xm) for all lensed galaxies in the 
sample. It is important to bear in mind that this fitting 
function for the SED is largely empirical, and the degree 
to which Td and (3 represent physical quantities is compli- 
cated by the spatial averaging over the entire galaxy and 
the degeneracy between a distribution in dust tempera- 
ture and a distribution of dust types (represented by ft). 
The formal errors for the fitted parameters (Td) should 
not be interpreted as errors on physical quantities, due 
to these caveats. 

Even though such SED fits could be obtained using just 
the photometric points, the addition of Z-Spec data not 
only strongly constrains the continuum slope, bu t also 
break s the degeneracy between Td and redshift (jBlainl 
I1999D , by independently determining the latter. Since Z- 
Spcc has determined the redshift, we are able to obtain 
Td from the continuum fit, which otherwise would con- 
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strain only the quantity T^j (1 + z) (e.g., lAmblard et al.l 
l2010h . This de generacy can lead to significant variations 
in the derived Td if the rcdshift is not measured indepen- 
dently. The implications of the continuum slope mea- 
sured by Z-Spec for the dust composition is left for a 
future work. 

Using Eq. [5] corrected for redshift and the derived Td, 
we can estimate the observed size of the dust emitting 
region. This solid angle will be affected by the lensing 
magnification factor. If the dust optical depth at submm 
wavelengths is low, as is often the c ase, fid will be corre - 
lated with r, and therefore with ft (jHughes et al.lll993l ). 
However, we break this degeneracy by fixing ft. The 
resulting fid ranges between 0.30 arcsec 2 for SDP.17b 
and 1.44 arcsec 2 for SDP.17a. These values may un- 
derestimate source sizes that are resolved in the SMA 
images with a resol ution of ~0.8 arcsec at 340 GHz 
(jNegrello et al.l l2010t) . Using the magnification factors 
from Table [3J the intrinsic size of the dust emitting re- 
gion will have an equivalent radius of 0.7 kpc for SDP.81 
and 1.3 kpc for SDP.130. Note however, that fid corre- 
sponds to the effective solid angle of the dust emitting 
region, such as the total area of small clumps spread over 
a larger region. In an image where these clumps are un- 
resolved, the total observed solid angle can appear to be 
larger. 

Having estimated the source size, the total dust mass 
follows from the relationship t{v) = K,(v)Md/ D\fld, 
where n{y) is the dust absorption coefficient k(v) = 
QA{v/2mGHzf cm 2 g" 1 (e.g., IWeifi et al.ll2007ah . and 
D^4 is the angular diameter distance. The dust mass can 
also be estimated in the optically thin limit (1 — e~ r ~ r) 
without the additional step of deriving fid , by substitut- 
ing t{v) directly in Equation [8] This is a good approx- 
imation at 250 GHz (1.2 mm), in the middle of the Z- 
Spec bandpass. Calculated in the optically thin limit, 
the dust mass is a robust estimate of the lower limit for 
the total dust mass in the galaxy, Mddim- Using the op- 
tically thin approximation and the 250 GHz flux density 
measured by Z-Spcc, we derive values for the magnified 
Mdjim of a fewx 10 9 M , as listed in Tabled If the dust 
is optically thick, as suggested by z;o=1300 GHz, the cal- 
culated Md.iim will underestimate the true dust mass for 
our galaxy sample by at most 30%. The dust mass is also 
inversely correlated with the assumed temperature, and 
will be underestimated when using the dust temperature 
corresponding to the peak of the SED. This temperature 
is likely too large to represent the bulk of the dust. As- 
suming that the 250 GHz flux is partially due to a dust 
component with a temperature as low as 20 K, and taking 
into account the optical depth corrections, we estimate 
that the total dust mass could be larger than Md,h m by 
up to a factor of ~4. To summarize, with good approx- 
imation, the true dust masses for these galaxies will be 
found in the interval [1,4] x Md,u m - The remaining un- 
certainties in Md.iim are mostly due to uncertainties in 
the expression for k,(v). Note that the quantity Md/fld 
is proportional to r and independent of temperature; for 
a given r, a lower limit for Md implies a lower limit for 
fid, but this limit for fid will decrease with increasing r. 

A more realistic approach is to fit the photomet- 
ric and continuum points with a library of SEDs, tak- 
ing into account the transmission curve of each instru- 



ment. We apply this method to our galaxy sample, using 
the SED libra ries of lCharv fc Elbazi ((200l (CE01), and 
iDale fc Heloul (120021) (DH02). The CE 01 templates have 
also been used bv lHwang et al.l (|2010| ) to fit both a low-z 
(z < 0.1) and a high-z (0.1 < z < 2.8) galaxy sample, 
of which the highest luminosity tail seems to have prop- 
erties overlapping the SMG population. Except for the 
subset of high-z galaxies with dust temperatures colder 
than ^90% of the local galaxies for a given luminosity, a 
subset that might be affected by blending, the CE01 tem- 
plate fits provide a good estimate for the total Lm in the 
high-z sample. For our lensed SMGs, we constrain well 
the peak of the SED and the dust temperature, and we 
are not in the regime where template mis match can have 
a big impact on the inferred Lm (see also iMurphv et al.l 

We find that the IR luminosities derived from the 
modified blackbody fitting are at most a factor of ^2 
lower than those when we use the SED libraries, and 
wi thin 20% from the L/b obtained assuming the models 
of Ida Cunha et al.l (120081) . cal ibrated for ULIRGs, with 
Av > 2 (jNegrello et al.l 20101 ). The variations between 
the values of Lm obtained by different methods reflect 
the systematic uncertainties in deriving this quantity. 
Similar underestimates have been found by others, and 
are due to the fact that the submm photometry does not 
measure the warm dust component of the SED, if one i s 
present (e.g., ISwinbank etlfl1l2T)10l : H vison et al.ll2010bl ). 
This is emphasized by the poor fit of the modified black- 
body curve to the Herschel-PACS data points, and the 
significant flux at shorter wavelengths predicted by the 
CE01 and DH02 models (Figure 0). Any derivation of 
Lm is model dependent, with the largest differences aris- 
ing from the presence of a warm dust component in the 
SED libraries. In Table [2] we list the Lm values derived 
from the SED template fitting method, as they represent 
a more accurate description of the total infrared energy 
output than the modified blackbody. The listed Lm arc 
calculated as the average between the values given by the 
best-fit CE01 and DH02 templates. We find that the two 
best-fit templates give values for the Lm within 15% of 
each other, both falling easily within our quoted error 
bars. On average, our Lm values are ab out 25% lower 
than those found bv lNegrello et al.l (|2010l) . but these dif- 
ferences are difficult to judge without data shortward of 
100 /im. Note that these values are rather smaller than 
typical IR luminosities of classical SMGs and unbiased 
sources detected by Herschel surveys. 

Except for SDP.17, we attribute all the submm flux 
density to the high redshift galaxy. The foreground lenses 
for the other galaxies in our sample have optical prop- 
erties consistent with being quiescent elliptical galaxies, 
and are therefore unlikely to have a significant submm 
emission. We have attempted a decomposition of the 
SDP.17 SED, using the two measured redshifts and a 
wavelength-independent scaling factor for each of the two 
components. The \ 2 value for the SED template fits 
is minimized when the observed flux density is split in 
half between the two components. This factor has been 
taken into account in Figure [9] and in deriving the Lm 
for SDP.17a and SDP.17b, as listed in Tabled However, 
the large dust mass inferred for SDP.17a could be an in- 
dication that the SED decomposition between SDP.17a 
and SDP.17b overestimates the contribution of SDP.17a. 



16 



Table 4 

Derived starburst properties and LTE parameters for the H-ATLAS galaxy sample. 



H-ATLAS 




^ ^ CO ,corr 




*SF (a) 


A r CO (b) 




/ S~2a 


TCO 


M CO {d) 


SDP ID 


(10 10 Kkms' 1 pc 2 ) 


(10 10 K km s" 1 pc 2 ) 


(10 11 M Q ) 


(10 7 yr) 


(10 17 cm" 2 ) 


(K) 


(io- 3 ) 




(10 6 M ) 


SDP.9 


13±3 


16±3 


2.1 


3.2 


23±4 


97±66 


0.3 


0.225 


1.6±0.3 


SDP.ll 


15±5 


18±6 


2.4 


2.1 


47±15 


48±8 


0.3 


0.929 


3.2±1.0 


SDP. 17a 


4±3 


5±3 


0.7 


11.1 


3±1 


160±90 


1.3 


0.005 


6.5±2.2 


SDP. 17b 


12±3 


16±4 


2.1 


3.7 


22±6 


80±17 


0.3 


0.114 


1.5±0.4 


SDP.81 


10±3 


15±5 


2.0 


2.1 


12±4 


62±8 


0.8 


0.133 


0.8±0.3 



Note. — The columns list: 1) the ID of the source in the SDP H-ATLAS catalogue; 2) the integrated brightness te mper ature of the lowest 
J CO transition measured, times the source area; 3) same as Column 2, corrected for dust absorption (see Section I4.2H ; 4) the molecular 
gas mass; 5) the gas depletion time; 6) the CO column density; 7) the CO excitation temperature under LTE; 8) the estimated beam filling 
fraction for the lowest J transition measured; 9) the optical depth for the lowest J transition measured in that source. The parameters in 
the last four columns have been derived in the LTE approximation; and 10) estimated total mass of CO gas. 

( a ) The errors for these parameters depend mostly on the uncertainties in the assumed conversion factors (see Section l4.2l l. 

( b ) The values displayed correspond to an intrinsic source diameter of ~2 kpc. The listed errors reflect the uncertainties in the measured 
integrated line fluxes. Other errors for this parameter, aside from the LTE model assumption, depend on our knowledge of the true source 
size. 

' c ' The formal errors bars underestimate the uncertainty in T ex , due to model assu mption s restricted to LTE. In non-LTE models, a large 
region of the parameter space is allowed, and T ex becomes J-dependent (see Section l4. 2.311 . 

' d ) Corresponds to the assumed source radius of 1 kpc, except for SDP.17a which would have an estimated radius of 5.5 kpc at z=0.942, 
estimated from the optical image. 



We estimate the star formation rates for our 
galaxy sample using the co nversion factor SFR(Mq 
yr-^)=1.5xlQ- 10 (L JR / L Q j ([Solomon fc Vanden Bout] 
120051) . similar to the iKennicuttl (|1998t) relation for 
a continuous starburst with a Salpctcr initial mass 
function (IMF). Since the selected galaxies are lensed 
by foreground object s with magnification factors ~10 
(jNegrello et al.l I2010D . the intrinsic IR and CO line 
luminosities will be ~10 times lower than the direct 
conversion from the measured fluxes. SDP.81 and 
SDP. 130 have magnification factors of 25 and 6, re- 
spectively, as derived from the best-fit lens model to 
the high-res olution sub-mm imag es available for these 
two objects (jNegrello et al.l l2010h . In Tables [2] and gj 
we left the quantities affected by gravitational lensing 
magnification unmodified, for reference, but the presence 
of this contribution is indicate d by the letter a in fr ont. 
Based on model predictions (jNegrello et al.l l2007t ). a 
typical amplification factor of 10 can be applied to these 
values. Once corrected for magnification, the infrared 
luminosities and corresponding SFRs are those typical 
of ULIRGs. 

4.2. CO Line Luminosities and Spectral Energy 
Distributions 

The measurements of CO lines reveal important in- 
formation about the physical properties and excitation 
conditions of the molecular g well as the total gas 
budget in these galaxies. These parameters can be used 
to investigate the link between star formation and gas 
properties. Higher gas temperatures and lower densi- 
ties would result in the increase of the Jeans mass, sug- 
gesting that star formation is bias e d towards high-mas s 
stars (e.g., lElmegreen et al.l 120081: [Klessen et all 120071 ). 
An increasing number of studies show that star forma- 
tion may proceed differently in merger/starburst systems 
versus quiescent /disk systems, the former being char- 
acterized by a top-heavy initial mass funct ion (IMF) 
(jWeidner et all 1201 It iHabergham et~al1l2010D . Such an 
IMF not only arises in dense starburst environments, but 



also has been inv oked to explain th e observed number 
counts at 850 /xm (jBaugh et alj|2005[) . A top-heavy IMF 
can arise in high-density material, shielded from far-UV 
radiation, but permeated by cosmic rays or X-rays, which 
heat the gas efficiently and generate cosmic-ray d om- 
inated regions (CDRs, IPapadopoulos et al.l I2010a|) or 
X-ray dominated reg ions (XDRs. Bradford et al.ll200l 
ISchleicher et~a l. 2010). The presence of a top-heavy IMF 
would have important consequences for the SFR inferred 
from the total Lir. However, the XDR signatures, suc h 
as highly-excited CO lines (e.g. Bradford et al.ll2009() . 
will likely indicate that the galaxy is dominated by the 
presence of an AGN, but not directly probe the IMF. 
Since the gas properties derived from the analysis of the 
CO lines are galaxy-averaged, only with multiple CO 
lines we can begin t o disentangle different PD R and XDR 
contributions (e.g. Ivan der Werf et all2tjTol) . which may 
help us characterize the star formation in these galaxies. 
Further understanding would require spatially resolving 
the star-forming regions, and probing the high-density 
star-forming gas with additional molecular tracers. 

In order to derive the physical characteristics of the gas 
in these galaxies, including the gas temperature, den- 
sity, pressure, and CO column density, we need mea- 
surements of multiple CO transitions, sampling the ro- 
tational ladder as fully as possible. The spectral line en- 
ergy distribution (SLED) for the CO molecule has been 
constructed in a few cases for near by and low redshift 
galaxies (e.g., iPanuzzo et al.l l2010h . In Figure ITUl we 
show the partial SLEDs for our galaxy sample, con- 
structed from the lines detected in the Z-Spec band- 
pass. This plot favors a distribution with the brightest 
lines between CO(5— 4) and (7—6), similar to the distri- 
bution observed for other SMGs and sta rburst galaxies 
(jWeifi et al.ll2007bt iDanielson et al.ll27LTol ). 

The shape of the line luminosity distribution does not 
reflect only the gas kinetic temperature, but also the gas 
density, and the effects of the optica l depth at the line fre- 
quency (jGoldsmith fe Langerl 119991 IPapadopoulos et all 
l2010bl) . In the optically thin limit, the CO column den- 
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sity scales with the absolute value of the line intensity, 
assuming that the source size and the magnification fac- 
tor are known. Under the assumption of local ther- 
modynamic equilibrium (LTE), all CO transitions have 
the same excitation temperature, T ex , also equal to the 
gas kinetic temperature Thini signifying that all rota- 
tional levels are populated according to the Maxwell- 
Boltzmann distribution at temperature T ex . In Sec- 
tion 14.2.21 we estimate these parameters by fitting the 
partial SLEDs, using the relationship between the inte- 
grated line brightness temperature, column density and 
excitation temperature, under LTE. Although this case 
is limiting due to the assumption of constant T ex for 
all levels, it is interesting to compare the predictions of 
this model to the more general non-LTE models, given 
its simple physical interpretation. In the non-LTE case, 
presented in Section 14.2.31 the models involve a larger 
number of par ameters, and are less wel l constrained. We 
use RADEX (|van der Tak et all [20071) to compute the 
brightness temperatures of the CO lines and estimate the 
likelihood distribution over the parameter space. These 
distributions allow us to assess if the available data are 
able to distinguish between the LTE and non-LTE mod- 
els. 

4.2.1. Gas Masses 

A useful quantity describing the CO lines is the 
velocity-integrated brightness temperature scaled by the 
area of the source, L co , in units of K km s _1 pc 2 . If 
the CO is thermalized and the lines are optically thick, 
L co will be the same for all rotational transitions for 
which the Rayleigh- Jeans approximation holds. In what 
follows, the brightness temperatures are computed in the 
Rayleigh- Jeans limit, and the values for L co are listed in 
Table E] both corrected and un-corrected for dust absorp- 
tion. Taking into account our estimate of the dust optical 
depth in Section 14. 1[ the observed brightness tempera- 
ture of the CO lines will be related to the intrinsic bright- 
ness temperature via the relation T obs = cxp(— Td)T mt , 
where = (v/v^Y . This correction will tend to boost 
the intensities of the higher J line s and drive the ex- 
citatio n of the gas higher (see also IPapadopoulos et aLl 
l2010bl) . The physical parameters of the gas derived in 
Sections 14.2 . 21 and 14.2 . 31 are also based on the absorption- 
corrected line intensities, and the effects of this correction 
are discussed as necessary. 

L co is traditionally derived from the CO(l— 0) tran- 
sition and related to the total molecular mass via 
the empirical relation M„ 



= aL co , where a is 



4.6 M q (K km s" 1 pc 2 ) -1 for the Galaxy ([Solomon et al.l 
fl997l). and 0.8 (K km s" 1 PC 2 )" 1 for ULIRGs 

dDownes fc Solomonlll998HTacconi et alJ l2008h. Follow- 
ing iSolomon fc Vanden Bout! (|2005l )~ we use the latter 



value for a and the L co for the lowest observed CO tran- 
sition to determine the gas masses (see Table [4]). This 
procedure assumes that all transitions from CO(l— 0) 
up to the lowest observed are thermalized, which might 
not necessarily be the cas e. A recent compa rison of the 
CO(3-2) and (1-0) lines (jHarris et al.ll2010l ) shows that 
the ratio of the brightness temperatures for these two 
lines averages to 0.6 rather than 1, due to the pres- 
ence of multi-phase CO gas. Moreover, the mid- J CO 



transitions do not account for the possible presence of 
a colder gas component, making the M gas derived in 
this manner a lower limit for the total gas mass in the 
galaxy. Assuming that the lines with J u >3 are ther- 
malized, corresponding to a warmer gas component, we 
apply this correction factor to the lowest CO transition 
measured, and obtain the gas masses listed in Table |H 
However, for subthermal excitation the ratio between 
the brightness temperatures of higher- J CO lines and 
CO(l— 0) could be even smaller. The value of ^c?o(i-o) 
for SDP.81 derived from the CO(l— 0) line intensity is 
1.8 x 10 10 K km s _1 pc 2 , after correcti ng for the lensing 
magnification factor (|Fraver et al.H20lTI ). which results in 
a brightness temperature ratio between the CO(7— 6) and 
CO (1—0) lines of 0.33±0.16. This also indicates that our 
conversion factors will globally underestimate the total 
gas mass. 

Using the SFRs derived from the IR luminosities, the 
gas reservoir probed by CO implies a gas depletion time 
[M gas / SFR] m these o bjects of ~ 10 7 years, similar 
to other known SMGs (jSolomon fc Vanden Boutl 120051: 
ICreve et al.l 121)051) . This can be interpreted as the star- 
burst lifetime under the assumptions of constant SFR 
and no gas inflow. Note that in the absence of differ- 
ential lensing, this estimate of the gas depletion time 
is independent of lensing magnification. We currently do 
not have enough data available to construct lensing mod- 
els and constrain the differential lensing for each of these 
sources. The star formation efficiency can be expressed 
directly in terms of Lm/ L co , without the need for a gas 
mass or SFR conversion factor. After accounting for the 
lensing magnification factor, Lm and L co for our sample 
follow the same relationship as other SMGs and ULIRGs 
(|Greve et al.ll2005l: IWang et al.ll2010l ). within the scatter. 

We derive an average molecular gas-to-dust ratio for 
the lensed galaxies of 127±50, subject to the caveats 
above: the gas mass is underestimated using the stan- 
dard conversion factor, and the dust mass is also 
underestimated by the single component model fit. 
The mean gas-to-dust ratio does not include the fore- 
ground SDP.17a, and is in agree ment with the val- 
ues found for other SMG samples dKovacs et al.l 120061 : 
iMichalowski et al.ll2010t iSantini et al.ll2010f ). Similar to 
the gas depletion lifetime, this ratio will be independent 
of magnification if we ignore differential lensing. 

4.2.2. LTE models 

The integrated line flux S^Av (in Jy km s _1 ) 
in the observer's frame is related to the velocity- 
integr ated Rayleigh- Jeans source brightness W(J) by 
(e.g., ISolomon et aLlll997t ) 



W(J) 



X 2 

A J,J- 



2 k ^^q. 



-S„Av 



(9) 



where VL S and VL a are the solid angles of the source and 
the antenna, respectively. W(J) is in units of K km s _1 . 
The last fraction represents the inverse of the beam fill- 
ing fraction. The contribution of the gravitational lens- 
ing magnification should cancel out in this expression, 
as it contributes to both S v and f2 s , but the true f2 s 
is not known. In principle, the same approach taken 
for the continuum (Section 14. 1[) could be used to de- 
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Figure 10. Spectral line energy distributions, uncorrected for 
gravitational lensing magnification. The Z-Spec measurements are 
shown connected by the bl ack histogram . The data point for the 
CO(l-O) line measured bv lFraver et~aH i20TH) in SDP.81 falls at 
the bottom of the panel, and is better seen in Figure [13] The red 
lines show the SLEDs predicted by the best-fit LTE model (con- 
tinuous), and the LTE models corresponding to the limits of the 
lcr standard confidence interval for T ex determined from the fit 
(dashed). The green line shows the SLED predicted by RADEX, 
using the parameters corresponding to the 4D maximum likelihood 
solution, as listed in Table 



termine the source size. However, such a fit requires 
a minimum of three parameters, and will not be well 
constrained by the number of CO lines in our SLEDs. 
In addition, the optical depth depends directly on the 
column density, and cannot be estimated independently, 
in the same way that the dust optical depth was de- 
termined by the continuum slope. We assume an in- 
trinsic source size of ^2 kpc. consistent with the an- 
gular dia meter of 0.2" of the SMG SMMJ2135-0102 at 
z=2.325 9 dSwinbank et~aT]|2010l: [Danielson et ail l2010t) 
used by INegrello et all (|2010[) for the SDP H- ATLAS 
sources, and similar to the size of the dust emitting re- 
gion found in Section|LlJ The corresponding beam filling 
fractions are listed in Table [4] As this source solid angle 
now represents the intrinsic size, and not the magnified 
one, we must correct the observed flux densities by the 
lensing magnification factors. We use the values listed 
in Tabic [2j when available, and assume a value of 10 in 
all other cases. The case of SDP.17a is treated differ- 
ently, as it is assumed to be a foreground galaxy, not 
affected by gravitational lensing. For the intrinsic size of 
SDP. 17a, we use a value of 1.54 ar csec 2 , which approxi - 
mates the size of the optical image. INegrello et al.l (|2010f ) 
identify two galaxies in the i-band image of SDP. 17 and 
fit both light distributions with the GALFIT software. 
As the presence of two galaxies could indicate a possible 
merger, we choose the source size of SDP. 17a to be the 
sum of the areas of these two galaxies. 

The distribution of the velocity-integrated bright- 
ness temperatures for the CO lines can be con- 
structed starting from the CO column density and gas 
te mperature, under the assum ption of LTE. Follow- 
ing iGoldsmith fc Langerl (|1999() , the velocity-integrated 
Rayleigh- Jeans source brightness is given by 



W(J) = N, 



hc z Aj.j-i 1 - 



where Tj t j-\ is the line center optical depth, and 
is the Einstein A coefficient for the transition. In LTE, 
the column density of molecules in the upper level, Nj, 
is related to the total column density N, by 



N, 



N 
~Z 



gje 



-Ej/kTe : 



(11) 



where Z is the partition function, Ej is the energy of 
level J, and g.; = 2J+1 is the degeneracy of level J. The 
line center optical depth can be expressed as a function 
of column density, temperature, and line width Av as 



Tj.J-l = A/,,7-1 



8tw 3 Av 



Nj{e 



hv/kT^ 



1). (12) 



Tj,J-l 



(10) 



We fit Eq. [TU]to the measured W{J) distribution, with 
the column density and gas temperature as free param- 
eters, and Av = 300 km s . We find that the best fit 
models have relatively low optical depths (< 1) such that 
the choice of the line width has only a small effect on the 
fitted parameters. For the lensed galaxies, the measured 
CO SLEDs and the range of SLEDs allowed by the for- 
mal lcr interval for the gas temperature are are shown 
in Figure [TUJ Since the CO lines are found to be close 
to optically thin in this model, the column density lcr 
interval would only scale these curves up and down, and 
not affect their overall shape. 

The SLEDs can be characterized by an overall scale 
and line ratios. The scale of the observed SLEDs is 
mainly a result of the CO column density and the beam 
filling fraction, while the line ratios depend on the CO 
temperature and gas (H2) density. The parameters 
in each pair are therefore largely degenerate and anti- 
correlated. This degeneracy is characteristic to CO and 
other molecular SLEDs, regardless of galaxy type. The 
last correlation (between temperature and gas density) 
only exists until LTE is reached, and the temperature be- 
comes fixed. By making assumptions on the beam filling 
fraction and gas density, we can place limits on the re- 
maining parameters. The error bars on the column densi- 
ties derived in this manner are correlated with the errors 
in the beam filling fraction, which are not known. Sim- 
ilarly, by making the assumption of LTE for all transi- 
tions up to CO (7-6), we are constraining the gas density 
to be greater than the critical density for this transition 
(n[H 2 ] > 3xl0 5 cm- 3 ). At densities n[H 2 ] > 10 6 cm- 3 , 
considerably larger than the average value observed in 
Galactic molecular clouds, all observed lines should be in 
LTE. Values of the gas density more typical for Galactic 
molecular clouds (10 3 -10 4 cm -3 ) correlate with higher 
gas temperatures, of a few hundred degrees, in order to 
reproduce the observed line ratios. 

The best-fit LTE CO column densities arc 
~fewxl0 18 cm" 2 , and the gas temperature ranges 
between 48 and 160 K, as listed in Table [H with the 
largest errors corresponding to the cases where only 
two CO lines have been measured. Note that these 
temperatures are derived after correcting the line fluxes 
for dust extinction, and are on average larger than the 
temperatures that would be obtained without correcting 
the line fluxes (between 41 and 115 K). However, due to 
the large errors in our measurements, these differences 
are not significant. 
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Taking into account the assumed source size, we es- 
timate total CO masses (Mco) of a fewxlO 6 M Q , or 
^10~ 4 of the total gas mass. This is consistent with the 
average relative abundance of CO, and would account 
for the entire molecular mass. However, since these LTE 
models imply very large pressures (~108 K cm ), the 
total molecular content of the galaxy would have to re- 
side in dense star-forming cores, which would account for 
the total CO emission. This suggests that the LTE pa- 
rameters cannot describe the overall average conditions 
of the gas in the galaxy. Other regions of the param- 
eter space are associated with non-LTE gas excitation, 
explored with the RADEX modeling in the next section. 

4.2.3. Non-LTE radiative transfer models of CO line 
excitation 

In general, the rotational levels of the CO molecule 
might not be populated according to a single tempera- 
ture, and the CO excitation temperature does not equal 
the gas kinetic temperature. By dropping the LTE as- 
sumption, we allow the excitation temperature to be a 
function of transition, being determined by the level pop- 
ulations for each line, while the kinetic temperature will 
be the global quantity describing the thermal energy of 
the gas. The level populations are found by solving the 
detailed balance equations including both radiative and 
collisional rates, and the output intensities are calcu- 
lated by solving the radiative transfer equations. Usually, 
these equations are strongly coupled, involving large spa- 
tial and frequency grids, and further complicated by the 
number of molecules and transitions involved. Simplify- 
ing assumptions are usually made to reduce the comput- 
ing time, depending on the problem at hand. 

We use RADEX to estimate the range of physical 
parameters consistent with the measured line strengths 
when dropping the LTE assumption. RADEX is a one di- 
mensional, non-LTE radiative transfer code, that solves 
for the level populations iteratively, employing the es- 
cape probability approxim ation for the radiative transfer 
(jvan der Tak et aLll2007t ). The medium is assumed ho- 
mogeneous and isothermal, and the number, type, and 
abundance of the participating molecules is selectable by 
the user. The input parameters are the kinetic temper- 
ature, Tkin, the number density of molecular hydrogen, 
npry, as the collisional partner, and the column densi- 
ties per unit line width of the participating molecules, 
only CO in our case. The background radiation field 
is the cosmic microwave background (CMB), redshifted 
according to the redshift of each galaxy. The output 
contains the predicted line excitation temperatures, op- 
tical depths, and line intensities. The output line fluxes 
are scaled by an additional factor 0, that represents frac- 
tional corrections to the size of the emitting region and to 
the gravitational lensing magnification factor. It would 
correspond to the area filling fraction of the emitting re- 
gion, if the size and lensing magnification factor of the 
source were known precisely. A value <f> > 1 would sug- 
gest that the assumed source size was underestimated. 
We compare the measured flux densities with the line 
intensities output by RADEX using the same values for 
source sizes, line widths, and lensing magnification fac- 
tors assumed in Section l4.2.2l for the LTE model. For the 
case 0=1 and n[H2]^ n cr i t for all transitions, RADEX 
will recover the LTE SLED as determined from T ex and 



N[CO] in Section |4"X21 as expected. 

We run RADEX for a range of input models, 
parametrized by Tki n , iV[CO]/dv, n[H2\ and 0, and com- 
pute the likelihood density fu nction for all model s follow- 
ing the method described in I Ward et al.l (j2003| ). Weak 
priors are set to rule out unphysical solutions, keeping the 
total molecular mass smaller than the dynamical mass, 
and the length of the C O column smaller than the phys- 
ical s ize of the galaxy (|Ward et al.l 120031 iPanuzzo et al.l 
120101 ). The dynamical mass cut-off is estimated choos- 
ing the line width of 300 km s _1 , and we require that 
the gas be self-gravitating (K vi r > 1, e.g., iScott et all 
l2011aHPapadopoulos et al.ll2l)l) Hi. We also impose a limit 
for the kinetic temperature at 3000 K, where collisional 
dissociation begins to rapidly destroy CO, weakly depen- 
dent on the gas density. 
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Figure 11. Contour plots of the (Tfefo, n[H2]) 2D marginal likeli- 
hood distributions, generated by a MCMC sampling of the param- 
eter space for RADEX models. The contours are in ncr-equivalcnt 
steps, enclosing 68.3%, 95.4%, 99.7%, and 99.99% of the probabil- 
ity, respectively. The dashed lines correspond to parameters that 
reproduce the LTE solution, and the dotted lines indicate the pa- 
rameters corresponding to the RADEX 4D maximum likelihood 
solution. Note that the 2D marginal distributions will not neces- 
sarily have the same maximum as the 4D distribution. The kinetic 
temperature is limited to < 3000 K, where collisional dissociation of 
CO becomes important. In the SDP.81 panel, the lighter contours 
show the probabilit y levels for a model including the CO(l — 0) from 
IFraver et al. (2011]). The parameters for this model are listed as 
model SDP.81* in Tablc[5] 

We map the surface of the likelihood distribution 
and determine the location of its maximum by running 
a Markov chain Monte Carlo (MCM C) algorithm, de- 
scribed in detail in IScott et al.l (j2011al ). Due to the large 
error bars and small number of data points, the afore- 
mentioned priors have only a weak effect on the final 
result, and mostly prevent the MCMC from spending 
time exploring unphysical regions of the parameter space. 
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Table 5 

Parameters used for the RADEX models shown in Figures j 101 and [T3l 



RADEX 




log(JVco) 


log(n[H 2 ]) 


4> 


log(P) 


log(M gas ) 


T 1 (d) 

CO total 


Model 


(K) 


(10 18 cm" 2 ) 


(cm" 3 ) 




(K cm- 3 ) 


(M ) 


(10 10 K km s" 1 pc 2 ) 


SDP.9< a ' 


99. 


18.79 


7.63 


0.44 


9.6 


8.96 


12.3 


68% credible region( b ) 


90.-1083. 


18.78-20.96 


5.39-7.97 


0.03-0.63 


7.44-10.39 


8.81-9.88 


5.5-37.3 


SDP.ll 


24. 


19.60 


7.35 


1.7 


8.7 


10.37 


15.5 


68% credible region 


24-612 


18.62-20.40 


4.20-7.93 


0.11-1.58 


6.63-9.93 


9.19-10.37 


0.9-32.0 


SDP.17b 


2833. 


19.73 


2.31 


0.91 


5.76 


10.24 


24.4 


68% credible region 


154-2884 


18.44-20.22 


2.31-6.38 


0.07-0.95 


5.76-8.69 


8.97-10.24 


1.9-38.6 


SDP.81 


375. 


20.09 


2.88 


0.26 


5.46 


10.04 


10.7 


68% credible region 


89-1416 


18.92-20.91 


2.88-6.13 


0.01-0.29 


5.45-8.28 


8.57-10.04 


0.2-31.2 


SDP.81*( C > 


453. 


20.07 


2.98 


0.18 


5.63 


9.85 


9.0 


68% credible region 


40-477 


19.66-20.35 


2.90-4.43 


0.17-0.74 


5.37-6.00 


9.70-10.18 


3.4-18.2 



Note. — The columns list: 1) the model notation; 2) the kinetic temperature T^i n (under LTE, 7fei„=T ex ); 3) the CO 
column density; 4) the density of H2; 5) <f> is an overall scaling factor, that would correspond to the area filling fraction 
if the intrinsic source size and gravitational lensing magnification factor were known exactly. This enters as the fourth 
unknown parameter in the maximum likelihood estimation; 6) the gas pressure; 7) the total gas mass in the beam; and 8) 

the Lju/ L C q total as a measure of the star formation efficiency predicted by each model, where Lqq is summed over all CO 
transitions in the model. 

( a ' These parameters correspond to the 4D maximum likelihood solution from an MCMC exploration of the parameter 
space with 10 5 iterations for each galaxy. Additional measured CO transitions would help rule out solutions with extreme 
temperatures and densities. 

( b ) This represents the smallest interval enclosing 68% of the marginal probability for each parameter. 

( c ) This second model for SDP.81 includes the CO(l— 0) measurement from lFraver et all ^2011f l. 

( d ) When derived from the integrated brightness temperature, a source radius of 1 kpc is assumed. 
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Figure 12. Same as FigureHTIfor the (N[CO], n[H 2 ]) 2D marginal 
likelihood distributions. 

The 2D marginal probability contours obtained from the 
MCMC algorithm are shown in Figures QT] and Q21 with 
the position of the 4D maximum likelihood indicated by 
the dotted line. Note that the 4D probability distribu- 
tions are highly non-gaussian, and therefore the coordi- 
nates of the maxima for the marginalized distributions 
in 2D do not match, in general, the paramters corre- 



sponding to the maximum of the 4D distribution. The 
set of parameters that maximizes the 4D likelihood for 
each galaxy is listed in Table [SJ and the line luminosities 
predicted by this model are shown in blue in Figure 1101 
The 68% credible regions are calculated as the smallest 
intervals containing 68% of the ID marginal probability 
for each parameter, around the value corresponding to 
the 4D maximum likelihood. In some cases, the credible 
regions for <fi suggest that the size of the emitting region 
could be larger than assumed for the LTE models, inter- 
preted as a larger area characterized by lower gas density 
and pressure than the LTE case. 

Due to the aforementioned degeneracies (see Sec- 
tion 14.2. 2[) , the product of the kinetic temperature and 
gas density on one hand, and CO column density and <f> 
on the other hand, are better constrained than individual 
parameters. These products are linearly proportional to 
the gas pressure and total gas mass, respectively, quan- 
tities which are listed in Table [5] 

The gas mass has also been derived in Section 14.2.11 
using two parameters: 1) the conversion factor a = 
0.8 between the gas mass and £^0(1-0) derived by 
iDownes fc Solomon! (|1998f ) from a non-LTE model, and 
2) the the scaling between the brightness temperatures 
of higher J CO lines and that of CO(1-0), R JA = 

Tp J ~ 1 /T 1 /., using the value R 3 1 = 0.6 ([Harris et al.l 
[2010f) . which we assumed to hold for higher J's. For com- 
parison, we can independently estimate both factors, a 
and Rj.i, using our best-fit non-LTE models. For a we 
get an average value of a — 0.46±0.24 for the whole sam- 
ple, not taking into account the error bars in the best-fit 
parameters. While for R31 we obtain an average value 
of i?3 ; i = 0.73, in relative agreemen t with the more rel i- 
able value of 0.64 ± 0.1 obtained bv lHarris etaLl (|2010f> . 
for the higher J ratios we have i?5 i = 0.75 (SDP.ll), 



21 




2 4 6 8 10 12 

Upper J level 



Figure 13. W( J) as a function of transition for four of the galax- 
ies in our sample. For clarity, the points corresponding to the 
same transition in different galaxies have been slightly offset left 
and right from the position of the exact upper J level. The triangle 
point re presents the intensit y of the CO(l— 0) line for SDP.81 mea- 
sured by Fravcr ct al. (2011). The W(J) distributions predicted by 
the LTF and non-LTE models are shown with a dashed and con- 
tinuous line, respectively. These lines emphasize the constraints on 
the allowed parameter space that can be gained by having measure- 
ments of both higher and lower J transitions. While the Z-Spcc 
data cannot clearly favor one of the models, the non-LTE model is 
superior when including the CO(l— 0) line for SDP.81. 

i? 6 ,i = 0.26 (SDP.17b), and for SDP.81 i? 7 ,i = 0.22 
and Rj i = 0.27 for the 2 non-LTE models, respectively 
(see Table [5]). This suggests that for higher J lines the 
Rj t i factor may be closer to a value of 0.3 for these 
excitation condit i ons. By comparison to the models of 
iNarayanan et al.l (|2009l ). values of 0.3 are marginally al- 
lowed, on the high end of the range. We emphasize that 
these estimates are strongly model-dependent, and only 
direct measurements of the CO(l— 0) lines would make it 
possible to both constrain the models and validate these 
values. It is becoming apparent that since most of our 
knowledge about local dust-enshrouded galaxies comes 
from the study of low-J CO lines, while at high red- 
shift the high-J CO lines are more readily accessible, 
we need to be able to measure both in order to make 
a direct comparison of the excitation conditions and gas 
properties. Important progress in this direction, by mea- 
suring the high-J lines in local gal axies, has been made 
with Herschel in rece nt years (e.g., IPanuzzo "eta!ll20T0l: 
IRangwala et al.ll201lD . 

The region of the parameter space that is most consis- 
tent with the observed line strengths is enclosed by the 
likelihood contours in Figures [TT] and [12] The likelihood 
space roughly splits into high density/low temperature, 
and low density/high temperature solutions. One addi- 
tional complication to the interpretation arises from the 
high dust optical depths, which lead to the suppression of 
CO lines with increasing frequency, and will cause an un- 
derestima te of the excitation tempera ture if unaccounted 
for (e.g., iPapadopoulos et al.ll2010'bl ). As mentioned in 
Section l4.2.2[ wc attempt to account for this effect by cor- 
recting the CO line strengths for dust absorption using 
the dust optical depths estimated in Section 14.11 How- 
ever, the likelihood distribution is relatively shallow over 
the whole region, reflecting the insufficient amount of in- 
formation in our data, and the likelihood contours are 



only marginally affected by this correction. 

To emphasize the insight gained by including addi- 
tional lines in the fit, we add to t he SLED of SDP.8 1 
the CO(1-0) integrated flux from iFraver et al.l (|2011l) . 
A likelihood analysis for the new set of lines results in the 
best fit parameters listed in Table [S] as model SDP.81*. 
The 2D marginalized likelihoods for this case are shown 
by the light grey contours in Figures [TT] and [12] The 
tightening of the likelihood contours is substantial with 
just one line added to the data, and the LTE region of 
the parameter space becomes less favored. However, the 
limitation of this model is that it assumes a single gas 
component, while most of the emission in the CO(l— 0) 
line could be originating from cold molecular gas. 

The constraints on the parameter space for the non- 
LTE models are weak, as expected given the limited sam- 
pling of the SLED and the large error bars, and cannot 
well distinguish between the LTE and non-LTE scenarios. 
The brightness temperatures predicted by both the LTE 
and non-LTE models are shown in Figure 1X31 emphasiz- 
ing the large deviations beween the predictions of the two 
models, especially for lower J transitions. The measured 
data points have been scaled by the lensing magnifica- 
tion factors listed in Table [5] when available, and by a 
factor of 10 in all other cases. This figure shows that 
the constraints on the model parameters can be tight- 
ened by measurements of lower J transitions, especially 
the CO(1-0) line. Even if most of the CO(1-0) emis- 
sion comes from a colder gas component, using this value 
as an upper limit will help rule out some regions of the 
parameter space, as in our example for SDP.81. 

The properties of the LTE and non-LTE models could 
be compared by calculating the total CO luminosity, 
summed over all transitions in the model, which is cor- 
related to the star formation rate and efficiency. Since 
the brightness temperature of the CO(l— 0) line tends to 
be lower in the LTE models, translating into a lower to- 
tal gas mass, the star formation efficiency, quantified by 
the Lm/L co ratio, will be higher in this case. As the 
temperature of the gas increases, more of the rotational 
CO lines become optically thick and high-J transitions 
start to dominate the gas cooling. Since the dominant 
cooling CO line is temperature-dependent, the total CO 
luminosity will be in general a better proxy for the total 
coolin g rate than t he lum inosity of a particular transi- 
tion. iBavet et al.l (|2009f) also find a strong correlation 

between the total L co and using a mixed sample 
of nearby and high redshift galaxies. If this relationship 
holds, we find that both LTE and non-LTE models are 
overpredicting the measured Lm, with a larger discrep- 
ancy in the non-LTE case. However, in both cases the to- 
tal L co integrated over all lines i s still consistent within 
the scatter with the IBavet et al.l (|2009D correlation, and 
therefore we cannot rule out cither scenario based on this 
comparison. 

Distinguishing between the different regions in param- 
eter space will clarify the state of the ISM in these 
galaxies, and thus their star formation histories. Specif- 
ically, hot/low-density gas may signal the action of a 
feedback process on star formation, increasing the Jeans 
mass. Other studies of high redshift SMGs find a 
warm CO component with n[H2] around 10 4 cm" 3 and 
temperatures between ~40 and 60 K (Riechers et al. 
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[20101 ICarilli et al.1 I20I0I: IDanielson et~al1 f2010h . a re- 
gion marginally allowed by our contours. However, a 
direct comparison with results obtained from SLEDs 
extending down to CO (2—1) becomes less warranted 
in view of increasing evid ence (e.g., iPanuzzo et al.1 
[20101 Bradford et "all T200l that the mid- and high-J 
CO lines are originating in some cases from a hot gas 
component. This high temperature/low density solu- 
tion has not been fully investigated, but recent stud- 
ies show that other CO SLEDs can be consistent with 
it (iScott etafl |2011at IPanuzzo et al.1 12010L iWeifi et al.l 
l2007aHAo et al.ll2008HBavet et alJ^OOaT The CO SLED 
in M82 is fit by a low-mass (~10% of the total) CO 
component with a kin etic temperature of almost 600 K 
(|Panuzzo et al.1 I2010I ). while solutions wi t h Ti- in of a 
fewx 100 K are found bv lScott et all (|2011a| ) ; iBavet et al.l 
(|2009f ). and ca n be allowed by the LVG models for IRAS 
F10212+4724 (lAq et al.l l200l and APM 08279+5255 
(jWeifj et al.ll2007al) . Similarly, the Herschel-SPIRE spec- 
trum of Arp 220 shows that the mid-J CO luminos- 
ity is dominated by a gas component with T~4350 K, 
which represents only ~10% of the total CO mass 
(|Rangwala et al.l 1*20 111 ). Such temperatures suggest en- 
ergy input from outflows or AGN activity. The presence 
of an AGN component in SDP.17b is supported by the 
relatively flat SLED from CO(6-5) to CO(8-7), similar 
to the Cloverleaf quasar o r Mrk231 (Bradford et al.ll200l 
Ivan der Werf e t al. 2010]), and the emission line of water, 
also observ ed in galaxies with an AGN co mponent, such 
as Mrk231 ([Gonzalez- Alfonso et al"1l2010h . 

5. CONCLUSIONS 

Far-IR / submillimeter-wave surveys are revealing 
submillimcter-bright galaxies from the first half of the 
history of the Universe by the tens of thousands, but their 
detailed study requires spectroscopic redshift measure- 
ments. We have studied a sample of the brightest sources 
and have demonstrated a new redshift-measurement 
technique with our broadband millimeter-wave grating 
spectrometer, Z-Spec. Z-Spec measures multiple rota- 
tional transitions of carbon monoxide, a major coolant 
of molecular gas in galaxies, and thus is not dependent 
on optical counterparts which are often absent or hard 
to identify, as is the case for these galaxies. We find red- 
shifts ranging roughly between 1 and 3, reaching back 
to an era when the Universe was 15% of its present age. 
Their fl uxes are proven to b e amplified by gravitational 
lensing (jNegrello et al.ll2fjTo| ) , making them ideal targets 
for spectroscopic follow-ups. From the observed CO line 
luminosities and integrated typical conversion fac- 
tors reveal that these galaxies each house roughly 10 10 
Mq of molecular gas, and have SFRs between 10 2 and 
10 M Q yr _1 , after correcting for lensing magnification. 
Regardless of the magnification details, we are clearly 
witnessing a rare episode of rapid star formation in these 
galaxies, since the timescale over which the observed lu- 
minosity can be generated by converting the inferred 
mass of gas into stars is only a few tens of millions of 
years (depending on the details of the star formation and 



the accretion of more gas), which is a small fraction of 
the Universe's age even at this early epoch. We estimate 
that the dust masses in our sample of lensed galaxies are 
around a fewx 10 s M Q , and the wavelengths correspond- 
ing to the peaks of their dust SEDs fall within a narrow 
range, between 73 and 92 /im in the rest frame. For this 
initial set of lensed submm galaxies both the dust proper- 
ties derived from the IR SED, and the physical conditions 
of the molecular gas probed by the CO lines, are broadly 
comparable to those in known SMGs dGreve et al.1120051 : 
iSolomon fc Vanden Boutl 120051 : iCasev et al.ll2009f ). with 
excitation temperatures in the 30-120 K range, and 
L co /L]r between 1 and 3xl0~ 3 K km s" 1 pc 2 /L Q , as 
measured from the mid-J CO lines. 

The partial SLEDs for the CO molecule constructed 
from the lines observed by Z-Spec cannot distinguish be- 
tween different models of CO excitation. The simplest 
assumption is that of local thermodynamic equilibrium 
(LTE), under which we can derive the gas column den- 
sity and excitation temperature. We find that the rel- 
ative line strengths can be reproduced by relatively low 
excitation temperatures (< 100 K), and optical depths 
(< 1). In the non-LTE case, other parts of the param- 
eter space are allowed, including higher optical depths, 
while measurements of the lower rotational transitions 
are essential in confirming such models. 

By being able to characterize galaxies that can be inac- 
cessible at other wavelengths, the combination of large- 
area submm surveys and spectroscopic follow-ups of the 
CO emission lines will lead to substantial progress in our 
understanding of high redshift galaxies and their evo- 
lution. These results suggest the possibility of a rapid 
growth in our understanding of high redshift star for- 
mation in highly dust-obscured galaxies, independent of 
identifying optical or radio counterparts, but enabled by 
strong gravitational lensing magnification. 
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APPENDIX 
A. 

Our redshift determination is based on denning the probability of false positives (or the false detection rate - FDR) 
in the absence of signal, and choosing the combination of estimators that produces the lowest FDR, and therefore the 
largest significance. In order to justify these choices and definitions, as well as the v~N normalization factor for E 2 , we 
will characterize and compare the distributions of the estimators, given that the signal Si in each channel is a normal 
random variable. This is the assumption of our simulations, which lead to the definition of the redshift significance. 
We verify that the distributions of the estimators are constant over the redshift range considered (0.5 — 6). 

A.l. Gaussianity 

All three distributions, denoted f(E 1 (z)), f(E 2 (z)), and f(E 3 (z)), respectively, are gaussian. This can be easily 
verified for f(E\) and f(E 3 ), since by definition they are constructed as linear combinations of normal random 
variables. E2 on the other hand, is defin ed, up to the normalization constant, as a sample median, which is a ce ntral 
order statistics. Numerous results (see !Shoracklll973t iRuvmgaart &: Van ZuijlerJ Il977t [Mason fc Shoracklll992l and 
references therein) show that order statistics, as well as the linear combinations of order statistics, of i.i.d. (independent 
and identically distributed) and non-i.i.d. variables are asymptotically normal. These conditions apply for the sample 
sizes n ~ 10 6 of our simulations, and therefore f(E 2 ) will also be well described by a gaussian distribution. 

A. 2. Expected values 

The expected values of the estimators should be at all redshifts for noise spectra, and should be largest at the 
correct redshift when lines are present in the spectrum. Taken independently or jointly, the values of these estimators 
determine the significance of the identified redshift. 

Let N(z) denote the number of CO lines falling in the Z-Spec bandpass at redshift z. In the noise simulations, for 
each channel i, 1 < i < N(z), the signal values Si are drawn from a normal distribution, with mean 0, and standard 
deviation equal to the noise value, cr^. Therefore, all Si/ Ui will be distributed as 7V(0, 1). 

It's easy to see that all our estimators have an expected value of in the absence of signal. In this case, E\ and E 3 
are just linear combinations of i.i.d. normal variables with mean 0. For simplicity, let's denote Si/o~i = Xi, where the 
Xj's are i.i.d. Af(0, 1), and re-write the definition of E 2 as 

E 2 {z) = y/N(z) x median(^), (Al) 

where A denotes the set 



A = {f ij \f ij = 0.5(x i + x j ),l<i,j<N(z),i<j}. (A2) 

The set A has M elements, with M = N(N — l)/2, and each element fcj = 0.5(xi + Xj) will be a Af(Q, 1/2) random 
variable. Since the expected value of the sample median is equal to the median of the underlying distribution, J\(0, 1/2), 
which is also (for a gaussian, the median is equal to the mean), the noise distribution of E 2 is also a gaussian with 
mean 0. 

If lines are present, let us assume that all Si (i.e. all channels containing a line) have the same mean So, and 
therefore are distributed as Af(So, erf). This is a simplifying assumption, which leads to a straightforward comparison 
of the estimators. In this case, we can also write the signal as Si = So + SSi, where SSi are Af(0, erf). In this case 
however, the Si/o~i ratios will no longer be identically distributed, each having a normal distribution with a different 
mean, Af(So/ai, 1). From the definitions given in Section [3~T1 we have for the expected values of the estimators: 



and 



where 



£(El (z)) = ^ = ^fL, (A3) 



£(E 3 (z)) = ^LJ2~ = ^0 <->, (A4) 

* / AT — * rr rr 



1 1^1 

(A5) 



° 2 >=ly° 2 

N ^ 



Note that the expected values of the estimators are calculated over all simulations, while the average of the CTj's is taken 
over the set of N(z) lines observed at redshift z. Applying the Cauchy-Schwarz inequality (< ab > 2 << a 2 >< b 2 >), 



24 



we have that < er 2 >>< a > 2 and V< a 2 > >< a >, and also 1 << a >< 1/cr > so that 1/ < a ><< 1/er >. This 
translates into 



S( El ( z) ) = -^£^= < ^ < V^VSo < - >= *(!%(*)). (A6) 

This proves that £3 has a larger expected value than E\, and will lead to a higher significance result than E\, when 
used independently. 

The expected value of £2 depends not just on some average of the noise values in all the channels where the lines 
fall, like E\ and £3, but on the distribution of these noise values. Depending on whether M is odd or even, either 
2, 3, or 4 channels will determine the value of median(A), and, depending on the noise distribution over the Z-Spec 
channels, the average noise in this subset of channels could be either larger or smaller than the average noise of all 
the channels containing lines. Due to this distribution of noise among the used channels, which is dependent on z, 
£(£2) is not consistently larger or smaller than £(£3) for any value of z, but we can set a limit for [£(£3) — £(£'2)1, 
independent of the values of Si and o~i, and therefore independent of z. 

By Chcbyshev's inequality, the distance between the mean and the median is always less than, or equal to the 
standard deviation, \mean(A) — median(A)\ < o~a- On the other hand, \J N(z)mean(A) is equal to £3: 

VN± £ 0.5(*< +Xj ) = ^^^(N - 1) E - = ^ E I (^) 

ij i i 

since each Xi will appear N — 1 times in the sum. The equality holds for any distribution of the signal-to-noise. As a 
consequence, the standard deviation of the sample mean for y N(z)A will also be equal to the standard deviation of 

£3. Intuitively, the normalization factor used for £2, \J N(z), allows us to apply Chebyshev's inequality in this way. 
For the standard deviation of A we have 

•i = 2(WTT) £ (« - ^ - to s2(i '» " W^) s « Al/n) - (A8) 

where s(xi) denotes the sample standard deviation, and x denotes the sample mean. Combining Equations IA7I and 
"T51 we obtain the constraint 

(E 3 (z) E 2 {z)f < ^"^ ^(IM), (A9) 

which shows that, since the sample A has a lower standard deviation than the original set {Si/ci} {{N — 1)/2(N+ 1) < 
1), the median of A will be closer than the median of {Si/o~i} to the average value of Si/oi. Therefore, £2 could be 
a better estimator than £3 for the case in which only a few of the N(z) channels are noisier than average (i.e. few 
outliers). 

The choice of any single estimator would be motivated by these individual properties. However, in the attempt to 
reduce the false detection rate ev en fu rther, we combine these estimators in pairs, by requiring that their maxima 
occur at the same redshift. Figure [ATI shows the combined FDR obtained for each of the three pairs of estimators, as 
a function of the estimator value (£1 for the (£i,£2) and {Ei,E 3 ) pairs, and £2 for the (£2, £3) pair). The FDR of 
a single estimator is also plotted with a long-dashed line, showing that for the same values of the estimator maxima, 
the FDR is lower in the combined case than in an individual case. The values of the Pearson correlation coefficients 
between these estimators, listed in the upper right corner, show that a lower correlation is associated with a lower 
FDR, since in this case the maxima are less likely to occur at the same redshift. The (£1, £2) pair is the one with the 
lowest correlation and FDR, and is the one used further in our algorithm. 

A. 3. Variance 

The last step in characterizing the properties of our estimators is deriving their variance. We will show that £1 and 
£3 have a variance equal to 1, while the normalization factor for £2 also brings is variance within a few precent of this 
value. 

Regardless of the expected value of the signal per channel, So, ^ar(£3(z)) = Var(Ei(z)) = 1, since 

£{El{z)) =1 + /^-, 

<cr > (A10) 

£(El(z))=l + NS 2 <-> 2 , 
a 

and Var(X) = £(X 2 ) — £(X) 2 , where £ (X) is given by Equations IA3I and IA4I for the two estimators, respectively. For 
the noise distribution, the variance of £3(2) also follows immediately from the fact that the variance of the sample 
mean of TV i.i.d. J\f(fio, c 2 ) random variables is Var(X) = o- 2 /N: 
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Figure Al. Comparison of the FDR's as a function of the pair 
of estimators selected. The long-dash line shows the FDR for a 
single estimator, E\. The notation Ek corresponds to E\ for the 
£2) and (Ei,Es) pairs, respectively E2 for the (E2, E3) pair. 
On the upper right corner we list the values of the Pearson corre- 
lation coefficient for the same pairs of estimators. Note that the 
derived FDR decreases as the correlation between estimators de- 
creases. The (Ei, E2) pair has the lowest correlation and also leads 
to the lowest FDR. 
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Figure A2. Upper panel: variance of th e med ian of set A, calcu- 
lated using the approximation in Equation lA18l with the correction 
for small N (solid line), and obtained by simulations (diamonds). 
For comparison, the dashed line shows the 1/7V curve. Lower panel: 
the standard deviation of E2, obtained by normalizing median(A) 
by the y/N factor. The solid line and the diamonds represent the 
analytic approximation, and the simulations, respectively, while the 
dotted lines are plotted to guide the eye. Note that the simulated 
points are in fact more linear than the semi-analytic formula, sup- 
porting an uniform y/N normalization factor, and the deviations 
from unity are only a few percent. 





Var(E 3 ) = —Var V — = NVar — V — = 1. (All 



It is trivial to see that for N = 2, the variance of E2 is also 1, since in this case E2 = E3. For larger values of N, 
the presence of correlations among the elements of A introduces important complications in deriving an analytic form. 
In fact, the covariance between any (fij,fik) = (0.5(xi + Xj),0.5(xi + x^)) pair is 1/4, and each fcj = (x{ + Xj)/2 
is correlated with 2N — 4 other variables. Therefore, the covariance matrix of A will have 1/2 on the diagonal, 
N(N -1)(N- 2) elements equal to 1/4, and the rest will be 0's. 

Up to the normalization factor, E2 is defined as the median of the set A, which is the central order statistic when 
M = N(N — l)/2 is odd, and a linear combination of order statistics when M is even (the average of the two values 
in the middle). An alytic expressions for the moments of order statistics have been derived for the case of i.i.d. normal 
random variables (|David fc Nagaralal 120031 : lTong||1990D. and generalized in a simple form only for the case of non- 
i.i.d. exchangeable normal random variables (|Owen &: Steckl 119621 : iTond fl990() . Exchangeable random variables are 
equicorrelated: the correlation matrix has all off-diagonal elements equal. However, the variables contained in the 
set A are non-i.i.d., and are not equicorrelated, except in the case when N — 3. For noise simulations, the fij's 
are drawn from the same distribution, but are correlated (ni.i.d), while if lines are present, the fij's will also have 
different means, becoming also non-idcntically distributed (ni.ni.d.). No analytic expressions for the moments exist 
for this cas e, and only a few genera l relations b etween the order statistics of such non-i.i.d. variables have been 
established dBalakrishnan et al.lll99l lTongH l990). We can show, however, that the variance of the median of A can 
be approximated analytically, and justify the choice of the normalizat ion factor for E2 . 

For a sample of i.i.d. normal random variables, it is well known (jCramerl 119461) that the sample median has an 
asymptotically normal distribution, with variance 

where Xm denotes the sample median, f(X) is the value of the distribution function at the position of the median, 
and M is the sample size. For a normal distribution N(uq,(Jq), f(X) = \ j ^/2-na} ) . Except for this asymptotic case, 
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there are no analytic forms for the moments of the sample median for the normal distribution, and the integrations 
have to be perform ed numerically. The value of Varad(XM) has been tabulated in the literature (jTeichroewl 119561 : 
ITietien et al1fl977h . We use the approximation 

, ~ . 4 2m m! 4 4 2m m! 4 , . , 

VarudiXM) ~ = — = 7^ TTTn i (A13) 

V ' (2m+l)! 2 27r/(A) 2 (2m+l)! 2 ' V ' 

based on the coefficient multiplying the exponential part of the distribution function, with m defined as M = 2m + 1 . 
Empirically, this expression offers a better approximation for small m's, than the one derived from the exponent of 
the exponential (given in Equation I A12[) .The last equality in Equation I A13I follows for an Af(0, 1) distribution. 

For a sampl e of non-i.i.d. equicor related random variables, with the same mean //o, same variance ctq, and same 
covariance Cn. lOwen fc Steckl (]1962[ ) showed that the variance of the median can be witten as 

Var nUd (X M ) « Co + (a 2 - C )Var ud (X M ) (A14) 

where VarudiX-M) is the variance for the sample median of M i.i.d. random variables distributed as 7V(0, 1). This 
relation provides an exact solution for the variance of E 2 in the case N = M = 3. Since the elements of A are 
Af(0, 1/2), and the covariance for any correlated pair is C A = 1/4, from Equation I A14I follows that 

Var A (X M=3 ) = i + (I - ^)var vld {X M ) ~ 0.361 « 1.08/JV, (A15) 

where we used Equation I A 1 3 1 as an approximation for small M. 

Wh ile an exact solution to the problem of non-equicorrelated variables is in principle possible jRawlings 197(1 iHilll 
I1976D , it involves a large number of integrations that ultimately have to be performed numerically. A similar situat ion 
arises in the study of genetic inheritance, and has led several authors (Meuwissen 1991: Ph ocas fc Colleaulll995D to 
develop approximate solutions, by assuming that all the variables arc equicorrelated, with a correlation coefficient equal 
to their ave rage correlation, and furt her refining this approximation using polynomial fits to Monte Carlo simulations. 

Following iPhocas fc CoTIeaul <fl995h , we can define an average covariance C e //, and assume that the sample A will 
behave as an equicorrelated sample with the new "effective" correlation. By definition, the variance of the sample 
mean is 

Var(X) = j^Yl Var ^ + E Cov(Xi,Xj). (A16) 

We define the average covariance, C e //, as the value that summed over all pairs produces the same total covariance. 
Since the covariance sum has M(M — l)/2 terms, we have 

TA .-, Mai 2 /Af(M-l)_ \ <7„ 2 M-1 ..„,. 

In order to derive C e //, let's remember that, since the mean of sample A is equal to E^/\fN (by Equation IA7|) . it 
will also have the same variance as E$/^fN, namely 1/N. By equating Equation IA17I with l/N, and taking into 
account that Ctq = 1/2, after some algebra we obtain C e // = 1/{N + 1). Substituting C e // for Co in Equation IA141 
the expression for the variance of the median for the sample A becomes 

1 N-l - 1 JV - 1 7T 

Var A (X M ) ~ h -; -Var iid (X M ) = h — ; 7-, r. (A18) 

Ay ' JV + 1 2(JV + 1) 1 ' N + l 2(N+l)(N + l)(N-2)' v ' 

where the last equality holds for large A^'s and can be approximated as 1/N + 0.57 /N 2 , and therefore has an overall 1/N 
behavior. For small JV's, we calculate the values of Var A {X.M) numerically, replacing Varud(XM) by Equation lA13l In 
this case, Var A {X.M) = 1/N + 0(1/N k ), where 0(1/N k ) has values less than 0.03. Thus, by semi-analytic arguments, 
Var A (X M ) a 1/JV. _ 

Even without additional polynomial corrections, the expression in Equation IA18I reproduces the actual variance 
within a few percent. We have checked this result by numerical simulations, obtained by drawing N(z) numbers 
from a Af(0, 1) distribution (corresponding to the Si/oi variables for the noise spectrum), constructing the set A, 
and taking its median. The expected va lue a nd variance of Var A {XM ) for each N(z) have been calculated from 10 6 
such samples. The upper panel Figure IA2I shows the analytic approximation with a solid line, and the simulated 
points as diamonds. T he 1/ N dependence is overplotted with a dashed line. For the analytic curve we have used the 
expression in Equation IA13I for N < 15 and Equation IA18I for larger N's. The bottom panel of the same figure shows 
a better comparison, obtained by multiplying the same curve and the points by N(z) and taking the square root. In 
this case, the values plotted represent the standard deviation of E 2 . From this figure it is apparent that, while our 
analytic approximation reproduces the behavior of Var(E2) within a few percent, based on the numerical simulations 
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Var(E2) is in fact even more linear than the approximation suggests, which further justifies the choice of the y/N as 
the normalization constant for E^. Even if the simulations do not asymptote exactly to 1, the difference is a constant 
multiplication factor, showing that there is no additional N dependence. 
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